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ABSTRACT 



The aim of this paper is to provide convenient predictor- 
corrector (P-C) methods for obtaining accurate numerical 
solution at a minimum cost to first order ordinary differen- 
tial equations (ODE). In pursuing this goal, a unified 
development of the most popular and efficient P-C methods is 
presented, which includes derivation of formulas and analysis 
of error propagation and numerical stability. Each method is 
then coded and programmed using the Fortran language. Com- 
parative analysis of the different P-C methods include both 
theoretical and numerical results. The numerical results 
were obtained by subjecting each method to a wide variety of 
test ODE, using a maximum of two corrector applications and 
a uniform series of step size values. By systematic compari- 
son of the performance of each P-C method the most convenient 
P-C ssts in terms of accuracy and minimum cost are established. 
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I. INTRODUCTION 



A linear first order ordinary differential equation (ODE) 
can be used as a mathematical model for a variety of phenom- 
ena, either physical or non-physical. Examples of such phe- 
nomena include the following: heat flow problems (thermo- 

dynamics), simple electrical circuits (electrical engineering), 
force problems (mechanics) , rate of bacterial growth (bio- 
logical science) , rate of decomposition of radioactive 
material (atomic physics) , crystallization rate of a 
chemical compound (chemistry) , and rate of population growth 
(statistics) . 

Most scientists, engineers or applied mathematicians, 
who have faced the problem of solving an ordinary differential 
equation numerically, are probably aware of the multitude 
of techniques available for such a problem. The abundant 
literature on the subject of numerical solution of ordinary 
differential equations is on the one hand, a result of the 
tremendous variety of actual systems in the physical and 
biological sciences and engineering disciplines that are 
described by ordinary differential equations and, on the 
other hand, a result of the fact that the subject is cur- 
rently active. 

The existence of a large number of methods, each having 
special advantages, has been a source of confusion as to 
what methods are best for certain classes of problems. It 
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is this observation which particularly motivated the writing 
of this paper; thus it is attempted to bring together the 
well-known predictor-corrector methods, which form a class 
of numerical integration methods for solving ordinary dif- 
ferential equations, in a consistent and comprehensive 
framework . 

This delimitation of the study to predictor-corrector 
methods is not without basis. Further research in the field 
showed that the predictor-corrector forms are the most 
efficient among the known integration methods in terms of 
speed and accuracy. Collectively, as a class of integration 
methods, the predictor-corrector sets are the best, but 
individually as a predictor-corrector set, the choice for 
the best method varies depending on the application. This 
paper atempts to compare the overall efficiency of all the 
well-known and most popular predictor-corrector sets by 
using a standard mode of application, the same series of step 
size values, and a set of test ODE's with many unusual and 
interesting features. To introduce the numerical experimen- 
tation of the different predictor-corrector methods, a 
comprehensive analysis of each P-C set starting with its 
derivation, error propagation, and stability is presented 
in detail. 

Numerical experiments are conducted using the IBM 360/67 
digital computer. The tremendous computational capability 
and speed of this computer offered an indispensable tool in 
conducting the experiment on a wide variety of test ODEs . 
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Single precision has seven digit accuracy and double pre- 
cision has fourteen digit accuracy available for computation. 

The paper starts with the description of the nature of 
ODE initial value problem. Then the various P-C methods for 
obtaining numerical solutions of the problem are enumerated. 
To lead into the derivation of formulas a brief review of 
backward difference operator and the well-known Newton Back- 
ward Formula is presented. In turn the different type of 
P-C methods are discussed followed by the analysis of error 
propagation and numerical stability. The stability bounds 
of the P-C methods were established through numerical 
experimentation. The next step is to subject each P-C 
method to a wide variety of test ODE. Then systematic com- 
parison of the performance of each P-C method is made. Con- 
clusions are derived and recommendations are made based on 
the analysis of numerical results obtained. Flow charts and 
computer programs are attached. 
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II. NATURE OF THE PROBLEM 



A. THE NUMERICAL PROBLEM AND NOMENCLATURE 

A linear first-order ordinary differential equation 
(ODE) of the form 

dy/dx = y'(x) = f(x,y) (1-1) 

with the initial condition 

y(x 0 ) - y 0 d-2) 

where f(x,y) indicates a differentiable function of the 
variables x and y is commonly referred to as an initial 
value problem. 

A typical elementary differential equations text presents 
several general classes of methods for solving a linear 
first-order ODE. The principal classes of methods are 
(1) variables-separable , or reduction thereto; (2) exact 
equations, or reduction thereto; and (3) solution by infinite 
series. The student is taught to apply the general method 
that appears best for the solution of the particular ODE. 

For example, the linear first-order differential equation 

dy/dx = xy (1-3) 

can easily be solved by the variables-separable method. 

This is accomplished by rewriting the equation in the form 

dy/y = xdx 

and integrating both sides to obtain 
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lny = x 2 /2 + c 



where c is an arbitrary constant of integration. The general 
solution can then be written as 



y = c,e 



x 2 /2 



Q 

where c^ = e . The general solution of such a linear first- 

order ODE consists of a family curves called the integral 

x 2 / 2 

curves. The family of integral curves y = c-^e ' that con- 
stitute the solution of equation (1-3) is illustrated in 
Fig. 1. 




Figure 1. Solution Curves of y' = xy. 



For each positive value of c^ a particular member of this 
family of curves is determined. 

A particular solution of equation (1-3) can be deter- 
mined if a condition on the solution curve is specified. 

For example, if it is required that the particular solution 
curve pass through the point (0,1), given by the initial 
condition 

y ( 0 ) = y 0 = 1, 
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is obtained since 



then the particular solution y 
C 1 = 1 * 




/ 2 



However, in real-life problems, many differential equa- 
tions that are encountered cannot be solved by elementary 
classical methods. In such instances, one must resort to 
numerical methods for obtaining one or more particular 
solutions of the initial -value problems. A multitude of 
techniques are available for solving such a problem numeri- 
cally. This paper considers only numerical methods using 
predictor-corrector equations. Even with this restriction, 
a large number of methods are in existence. Literature 
research in this subject, however, reveals that the most 
popular and efficient methods are the following: 



1. Euler Predictor-Corrector Method 

2. Milne Predictor-Corrector Method 

3. Nystrom Midpoint Predictor-Euler Corrector Method 

4. Hermite Predictor-Milne Corrector Method 

5. Milne Predictor-Hamming Corrector Method 

6. Adams Predictor-Corrector Methods 



These methods, each having special advantages, have been a 
source of confusion as to what methods are best for certain 
classes of problems. This paper will attempt to present a 
comparative analysis of these various predictor-corrector 
methods . 

In attempting to compare the various predictor-corrector 
methods, such questions as the importance of the number of 
function evaluations, the step size to be used in a numerical 
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calculation, the computing time required for each method, 
the influence of truncation and round-off error and the 
stability of various predictor-corrector modes of computa- 
tion will be considered. 
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Ill . DERIVATION OF PREDICTOR-CORRECTOR EQUATIONS 



A. BACKWARD DIFFERENCE OPERATOR 

To lead into the development of certain equations of 
major importance in numerically solving ODE, a linear back- 
ward difference operator V is defined by the following 
equations 



In general, these formulas start with a given sequence of 
the y (x ) = y^ the ordinates at equally spaced x n = x 0 + n ^> 
h is a constant. 

B. NEWTON BACKWARD DIFFERENCE FORMULA 

Using these operators the well known Newton backward 
difference interpolation formula is: 




( 2 - 1 ) 




3 2 2 

V y = V y - V y , 
7 n 7 n 7 n- 1 






( 2 - 2 ) 



+ 



q(q+l) 



ISLtH:- 1 -! y n y 



n 



where q = — — 
h 
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This formula is obtained by fitting a polynomial to the set 
of points (x- ,y ) at equally spaced (x n >. (Actually 

the formula is a polyn mial in a because of the change of 
variable.) Fitting is taken to mean that the polynomial 
passes through the points ( x n >y n )* Since there are 
(n+1) points (y ,y ^ , . . . , y ) , this means that if the solution 
function y(x) is a polynomial of degree n or less, the 
formula above will be exact. When y(x) is not a polynomial, 
or is a polynomial of degree greater than n, the formula 
will not be exact (except at the points themselves). In 
other words an error will be made because a finite rather 
than an infinite series is used. In general form this error 
will appear as 

T = C h n+1 y tn+1] (5) 
a a . 1 . 

where C is a function of a and y n+ ^ (?) is evaluated at 
some ?, x <?<x. T is commonly referred to as the 
truncation error. 

C. INTEGRATION FORMULAS FOR ODE 

Integrating both sides or equation (1-1) between ( x n >y n ) 

and (x .,,y„.,) there results 
n+ 1 ’ ' n+1 J 

x n+l 

y n+l = y n + / f(x,y)dx (2-3) 

x n 

X . i 

n+ 1 

= y n + / y ’ (x)dx 
x n 

1 

= y + h / y 1 (a) da 
n 0 
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Attaching the index sequence of n = 0,1,2,... to this equa- 
tion and noting that y Q is known at x , ( 1 - 2 ), it becomes 
apparent that (2-3) can be used recursively to generate 
yi ,y 2 ,y 3 > • • • as long as there is some way to evaluate the 
integral. In other words the solution of the ODE is con- 
verted to evaluating an integral. 

Derivations of many of the equations used to represent 
(2-3) now follows: 

a) Finite-difference table of backward differences. 
Assuming that (x^y^), (x£ ,y 2 ) • • • (x r , y n ) are given, and 
computing y' = f(x^,y^) for i = 0 , 1 ,... n the following table 
is constructed. 



x 



1 



X 



2 



Vyj 

V 2 y^ 

vy 2 



V q y ' 

7 n 



n-1 



y n-l 



n 



y f 

J n 



* y n 



V 2 y ’ 
7 n 



These values will now be used to continue the solution in 

t h 

computing y In general, by retaining q^ differences 
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in the Newton Backward Formula (NBF) (2-2), the following 
finite series results. 



y' = y' + aVy 1 + V 2 y* + ... 

1 n+a 'n 'n 2 ! n 

+ a(a + l) . .a(a+q-l) y q , 
a ! y n 



(2-4) 



Using this interpolating polynomial to approximate the inte- 
grand in (2-3) [replace y' (x) by y* ] yields: (complete 

1 1 ’ 

details of integration are presented in McCalla [Ref. 1]) 



h i 

y ., = y + h E a.V y , with a = 1 
'n +1 7 n i 'a’ o 



where 



u = f 1 da> 

1 n i ! 



for i > 0 



(2-5) 



Equation (2-4) is in extrapolation form as can be seen from 
the table of differences, where the end point x n+ i is 
excluded from the interpolating points. 



The error term associated with truncating after the q 



th 



V is 



t = h q * 2 / 1 SLC«±n.-.^(i»i a l y [q* 2 ) (5) dc 
a „ (q + l) I ’ 4 

r + ? 1 

But since the coefficient of y q does not change in sign 
in [o.lj, it is possible to write 



T = a , h q+ 2 y tq+2] ( 5 ) 
a q+l J J 



(2-6) 



By actually calculating the a^ in (2-5) one arrives at 
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y 



n+ 1 



= y 



n 



h[l + f 



V + 



5 

IT 



V 2 + 



3 3 

_ v 
8 



251 

72U 



V 4 + 



Jy' 

7 n 



Note that by truncating this series at the first difference 
then (for q = 0) , 

y i = y + hy' with T = \ h 2 y [2] (£) • (2-7) 

7 n+l 7 n 7 n a 2 7 v 7 

which is the well known Euler's formula. Adopting the con- 
vention that 

T a = T(x,h) 

since a is a function of x and h. Also 



T(x,h) = | h 2 y t2, (S) 



can be written as 

T(x,h) = 0 (h 2 ) 

2 

indicating that the truncation error is of order h . 

Actually (2-3) can be generalized by rewriting it in the 

form 



n+1 



= y +h/ y' da 

'-n + rv 



n+a 



( 2 - 8 ) 



where r is any positive integer. For example, the case r =0 
is the one already discussed. Following the same procedure 
as previously presented, the results for r = 1,3,5 are 
obtained by substituting (2-4) into (2-8) where 



a 



i 



/ a(a + l) . • . (a+i-1) da ^ 
-r 1 * 



for i 



1,2 



q- 



Actually calculating a^ for different values of r there 
results 
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r = 1; y n+l = y n-l + h[2+0V+ | v2+ | v3+ |f v4 + ( 2 - 9 ) 



]y ' 



n 



r = 3 



y 



n+l 



= y n-; 



h[4-4V+ § 



2 3 

V +0V + 



ii 

45 



V 



( 2 - 10 ) 



]y 



I 

n 



r = 5; y xl = y _ + h [6-12V+15V 2 -9V 3 + ~ V 4 + 0V 5 + (2-11) 
/ n+l / n-5 10 



The error term associated with truncation of these formulas 
is more complicated to calculate because the integrand 
changes sign in [-r,l]. For special cases however, as in 
the above for r = 1,3,5, the coefficient of the r 1 "* 1 dif- 
ference is always zero. Thus it is only necessary to use 
the r-1 differences of y' in order to calculate a result 
whose truncation error is of the order of r+2 in h. For 
r = q = 1 

>Vi - Vi ♦ 2h K C 2 ' 12 ’ 



with 



T(x,h) 





(O • 



Equation (2-12) is known as the Nystrom Midpoint formula. 
Similarly for: 



r=3; q= 3 ; y 



n+l 



= y n-: 



4h [y ' -Vy • 
; n 1 n 



2 2 
- V v ' 1 
3 y n J 



(2-13) 



with 
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T(x,h) = || h 5 y [5] CO 



r = 5 ; q = 5 ; 



y = y r + 6h[y' - 
'n+1 'n-5 7 n 



2Vy ' + — V 2 y' V 3 y' + — V 4 y'] 
y n 2 y n 2 y n 20 y n J 



with 



T(x,h) = ^ h 7 y [7] CO 

Generally, formulas derived using the extrapolation form of 
the NBF are referred to as open, explicit or predictor 
equation because y n+ -^ occurs only on the left hand side of 
the equation. In other words, y n+ ^ can be calculated 
directly from the right-hand side values. 

b) Extending the table of backward differences by one 
line to include x n+1 as an interpolating point yields 



y' 

' o 



y i 



Vy' 



2., i 



v y 



V q y ' 
1 n 



V 7 n 



c y ' 

n 1 n 



V y A+l 



c n+l y n+l 



v 2 y 



n+1 



V 4 v* 
v y n+l 



NBF now in interpolation form yields a finite series as 
follows : 
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(2-14) 



y 



n+a 



= y 



n+ 1 



+ (a- 1 ) y 



n+1 



(a-1) (a) 

21 



"y * 

' n 



n+1 



(a-1) (a) (a+1) . . . (a+q-2) q , 

q* Y n+1 

Repeating the procedure used in the extrapolation form 
(2-4), analagous results for r = 0,1,3 and 5 are 



r = 0 ; 



y 



n+1 



= y 



n 



h[l 
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720 



+ 



• • • ] y« 



n+1 



(2-15) 



r = 1; y 



n+1 'n-1 



y n l + h [ 2 - 2V + i V 2 + OV 3 



-1- V 4 
90 



• ]y: 



n+1 



(2-16) 



r = 3; y +n =y + h [ 4 - sv + ^ v 2 - 4 v 3 + ii v 4 - ov 5 
’ 'n+1 'n-3 33 



45 



r « 5; y - = y , + h[6 - 18V + 27V 2 - 24V 3 + V 4 

’ 'n+1 'n-5 10 

- to v5 + (2 - 185 

Some interesting and useful results may be obtained by trun- 
cating after a certain number of differences. In the r = 0 
case truncation after the first difference (q=l) yields 

y n+l = y n * 7 ty ibl * y A ! (2 ' 19) 

with 

T(x,h) = ^ y [3] (O 
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which is called the modified Euler formula. For r = 1 , 
truncating after the second difference (note that the third 
difference is zero) yields 



>Vi = y n -i + - 



h 

3 






, + 4 y' + y 
n+1 ^n y 



A-] 



( 2 - 20 ) 



with 



T(x,h) = ^ y t5! ( 5 ) 



which is popularly known as Milne's corrector. Formulas 

resulting from the use of NBF interpolation form are referred 

to as closed, implicit or corrector equations since y n+ ^ 

occurs on both sides of the equation. In other words 

unknown y cannot be calculated directly since it is con- 

tained within y' 

' n+ 1 

c) In sections a and b the equations were obtained 
directly from manipulating interpolation or extrapolation 
polynomials using NBF. However, another method, called the 
method of undetermined coefficients, applied to the gen- 
eralized linear, K-step dif ferential -dif ference equation 
with constant coefficients 



n+1 



K 

= Z 

i=0 



K 



cx • y 
i/n-i 



+ h 



Z 3-y' . 

i-1 1 n_1 



( 2 - 21 ) 



will yield all the formulas derived so far and at the same 
time can be used to obtain all other predictor- corrector 
equations by suitable choice of the parameters and 3^. 
But since the previous formulas derived are sufficient for 
applications and only the Hermite predictor is of interest 



25 



using the present method, the tedious and long process of 
derivations, which can be found in most numerical analysis 
books [Refs. 1,2,3], are left out and only the results for 
the Hermite predictor from [Ref. 4] are presented. There 
it was found that 

a = -4 a n = 5 a-, = a, = 0 

o I Lb 

*0 = 4 ^1 = 2 ^-1 ^2 ^8 = 0 
yielding 

y ,, = -4y + 5y , + h [4y' + 2y' ,] (2-22) 

7 n+l 7 n 7 n-l 7 n 7 n-l v J 

with 

T(x,h) = ^ y [4] (?) 

d) It is of interest to note that the numerical methods 
developed for solving the initial value problem differ in 
their requirements of the number of starting values. For 
example, Euler's formula (2-7) needs only the initial con- 
dition y (x q ) = y Q and h to start the solution. Methods that 
determine y n+ ^, when only one point ( x n >y n ) and step size 
h are known, are commonly called one-step method. 

On the other hand, the Nystrom Midpoint formula (2-12) 
needs starting values ( x n _i>y n -i) and ^n’^n-^ to contanue 
the solution. Methods that require step size h and more 
than one point (x ,y ) , (x ,y ,).... in order to compute 
y n+ l are called multi-step methods. 

It is clear then that a multi-step method requires a 
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one-step method to provide the necessary starting values 
in its computation. 



Literature research showed that the Runge-Kutta method 
is the best among the available one-step methods. Hence the 
fourth order Runge-Kutta method based on the Taylor's series 
expansion truncated after terms of the fourth order (deri- 
vation of this method is given in Ref. 3) was used as a 
starting method for the multi-step methods considered. This 
has the form 

K 1 +2K 2 +2K 3 +K 4 

>Vl = >^n + 6 

where 

K 1 * 

K 2 " hf ( x n + y n * 

K 3 = hftx n * I’ y n + I 1 ’ 

K 4 = hf(x n * h, y n * X 3 ) 



to solve the initial value problem 



y* = f(x,y) with y(x Q ) = y Q 



yielding as many starting values as needed. 
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IV. PREDICTOR- CORRECTOR METHODS 



A. DEFINITION 

Algorithms that use both an explicit formula and an 
implicit formula are called predictor-corrector methods. 

The solution approximation computed by the explicit formula 
is denoted y^ + ^ and is called a predictor. This predictor 
is then used initially in the right side of an implicit 
formula that computes a corrector Y n+ ^* The implicit 
formula can be repeatedly applied, using y^ + ^ from the 
preceding iteration on the right side and computing a new 
y 0 ^, on the left. 

To illustrate the procedure, the Euler P-C algorithm 
(one-step method predictor) and the Nystrom P-C algorithm 
(multi-step method) are used. 

Given the ODE 

y' = f(x,y), with initial solution point ( x n >y n ) 

1) Euler P-C method uses the equation (2-7) as predic- 
tor, and (2-19) as corrector 

■ >Vi = + hy A f2 - 23) 

>vi ■ * 7 (2 ' 24) 

Since (2-23) requires only one solution point ( x n ,y n ) and 
step size h (assumed to be constant for specific application) 
which can be chosen arbitrarily (the choice actually affects 
the stability and error propagation of the method as will 
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be illustrated later), y' = f(x ,y ) can be evaluated. Hence 
(2-23) can be computed yielding the result y^ + ^ (p means 
predicted value) . Then the derivative of yP + ^ i- s computed 
yielding the result y^ + ^: 

y ' , = f (x , ,y^ , ) . 

7 n+l v n+l ,7 n+l 7 

Next (2-24) is used to compute a corrected value called 
y^ + l (c means corrected value) 

y n + l “ y n + 7 [ y ^l + ' 

At this point there are two alternatives. First, the compu- 
tation can be terminated using Y n+ ^ as the true approxima- 
tion value for y Then the computation continues to 

determine the next solution point y n+ 2 by repeating the 
entire procedure. Second, the corrected value, y n+ ^ of 
y n+ ^ may not be acceptable in the sense that the equality 
of both sides of (2-24) is not satisfied to as many digits 
as may be desired. In other words, the corrector has not 
converged to the desired accuracy (convergence of the 
corrector formula is an important topic and discussion of 
this aspect is deferred to a later section) . The next step 
then is to improve the value y n+ -^ on the left side of (2-24) 
by taking the derivative of y n+ ^> calling the result y^ + ^, 
and then calculating an improved value of y n+ ^- This is 
repeated until convergence occurs to as many digits as 
desired. This convergence term, designated as EPS, is 
tested against the absolute difference between the last two 
approximations . 
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2) The Nystrom P-C method has (2-12) as predictor and 
(2-19) as corrector. 



It is clear that before (2-25) can be evaluated another solu- 
tion point (x ,,y _,) is needed in addition to the given 



solution point (x ,y _,) of the ODE. Another one step method 



is therefore needed to provide the additional starting 
values. The Runge-Kutta method can be used to accomplish 
this. Once these starting values are known, the same pro- 
cedure as in Euler P-C method is followed to compute the 
next point using equation (2-25) as predictor and (2-26) as 
corrector . 

B. A SIMPLE PREDICTOR- CORRECTOR SET 

The predictor- corrector methods discussed so far do not 

take into account the truncation error incurred in retaining 
t h 

r differences of the infinite series of the NBF. The 
application is straightforward, predict then correct, and 
as such it can be called a simple predictor-corrector set. 

C. MODIFIED PREDICTOR-CORRECTOR SET 

The P-C set given by (2-12) and (2-19) has a local trun- 
cation error of the same order as the set (2-23), (2-24) 
(T(x,h) = 0(h^)). Hamming [Ref. 4] has used this feature to 
expand the method of calculation. The local truncation 
errors for the predictor (2-12) and corrector (2-19) are 




(2-25) 




(2-26) 
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(2-27) 



where 



T < x - h )p - Vi.p * f>' 131 « P ) 



x , < £ < x ... 

n-1 p n+1 



and 



T C*- h ’c = Vl, 



■ iy r Isl «c> 



(2-28) 



where 



x < £ < x 

n n+1 



It follows that the exact value y(x n+ ^), 0 f y at x j, is 
given by 



>’( x „ + l) - vl 



♦ — y [3] (C 1 
n+1 3 y *^p J 



y(x n *i) ■ y 



n+1 



■^y [3] (c ) 

12 7 K c J 



From these equations, one obtains 



y C - lL_ y[31 (P ) = yP 
y n+l 12 7 u c J y 



+ iLl y^ ( E ) 
n+1 3 7 U P J 



or 



c p 

y n+l " y n+l 



r yI3]( V + ti yt3] (c c ) 



(2-29) 



Now assuming that 



y t3] (C p ) « y t3] (? c ) = y [3] (O 



Over the interval of interest (i.e., the third derivative 
does not vary greatly over the interval), (2-29) becomes 

yS*i - ■ n 1,3 y [I, «> (2 ' 30) 
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where 



Vi < 5 < Vi 



Even at this point in the development some interesting fea- 
tures can be identified. First one notices that (2-30) can 
be obtained only in case the order of the predictor equation 
is equal to the order of the corrector equation. Second, 
there is now a measure of the local truncation error in terms 
of yC +1 an d yj +1 which are explicitly calculated in the P-C 
algorithm. This may be seen more explicitly by comparing 
(2-27) with (2-30). 



c p 

>vi - y 



5 ,3 3 r >.n 

n+l 12 h y (5) 



5 r l u 3 3 r > r ^ 

= j (3 h y (£)) 



5 ^ 

T n+l,p 



Similarly comparing (2-28) with (2-30) 



c p 

>Vi - Ci 



-5 T 



n+1 ,c 



(2-31) 



(2-32) 



Since the local truncation error can now be estimated, the 
next question is how to use this information to improve the 
predicted and corrected values. Considering the predicted 
value first, it can be seen from (2-31) that 

T = — fy c - yP ) 

n+l,p 5 l/ n+l 'n+l' 

Assuming that the local truncation error remains approxi- 
mately constant over two steps, then 
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T 

n+l,p 







The predicted value can then be improved or modified by 
adding the term (y^ - y p ) to y p +;L using only information 
calculated previously. 

The corrector can be modified in essentially the same 
manner since 



X = - i fy C 

n+l,c. 5 



- yP ") 
n+1 x n+l-' 



can be added to y n+ ^ to improve this value. On this basis 
the overall P-C calculation can be written in the following 
algorithmic steps: 



Predict : 


P n+1 


= y n- 1 


+ 


2hy; 




Modify : 


m n+l 


= p n+l 


- 


I ^ p n' c n 


) 


Reevaluate 
Derivative : 


m n+l 


= f(x J . 1 
^ n+1 


’ m n+l^ 




Correct : 


c n+l 


= y n + 


h 

2 


(m' + 

^ n+1 


fn) 


Modify : 




c n+l 


+ 


I ^ P n+1 


' c n+l> 



Then iterate the corrector to convergence. To simplify the 
presentation, the symbols m , p and c have been used. 

The use of the modified predictor-corrector set seems 
to be very attractive but it must be noted that this method 
has two limitations. First, it rests on the fact that the 
original P-C equations have equal-order truncation errors, 
and second, it neglects round-off error. 
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D. HAMMING MODIFIED PREDI CTOR- CORRECTOR SET 



Hamming used equations (2-10) and (2-16), which are known 
as the Milne P-C set when both are truncated after the third 
differences, and he revised the corrector to remove the severe 
stability problems, which will be seen later. 

To develop the revised corrector, Hamming started with 
the generalized equation C2 -21) and obtained 

>Vl = <Vn * “l y n-l * Vn-2 + ht6 -l>'A+l + 6 o y A 

* 6 l>n-l ] • 



Using the method of undetermined coefficients [complete deri- 
vations can be found in Ref. 4] , he obtained 



°o = I ( 9 - 9 “P 

a 2 = ’ "8 ^ 1_a l^ 



6 0 - £ C9 * Tap 

6 1 “ h ( ' 9 * 17a P 



with 



T(x,h) 



(-9 + 5a x ) 
360 



h S y E5] 



(5) 



with as a free parameter. Hamming then considered the 
values of a 0 ,° l 2*'* , ^l as f unct i° ns of a i for a ± ~ 1» 9/17, 
1/9, 0, -1/7, -9/31, and -6/10. After looking at the result- 
ing coefficients in terms of equal magnitude, number of 
zeros, etc., the case = 0 was chosen as the best value 
since it yields a greater region of real negative stability 
and a wide range of relative stability. It can be verified 
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that the case = 1 yields the Milne^ corrector (2-10) 
truncated after the third difference. With a-^ = 0 the 
Hamming predictor-corrector set is of the form 



Cl - y n -3 * I h (2 ' 33) 

which is (2-10) in the case r - 3, q = 3, of the NBF extra- 
polation form, with 



and 



T(x,h) = if h 5 y l5] (£) 

Cl - f ‘%-yn-z 1 * I h iy; + i * (2- 34 ^ 



with 

T(x,h) = ^ y [5] (?) . 

Next, using the procedure discussed in the previous section 
on the modified predictor-corrector set, Hamming develops 
his modified predictor- corrector algorithm as follows: 



Predict : 


Pn + 1 = 


Vs * 


3 h ^ 2y n" y n-l +2y n-2^ f 2 ' 35 ’ 


Modify: 


m n+l 


P n+ 1 ' 


'inHPn-S^ 


Reevaluate 
Derivative : 


m n+l 


ffx n*l 


,m n+P 


Correct: 


C n+1 = 


F '% 


-y -,] + h [m' +2y'-y' ,] 

7 n-2 8 n+1 7 n 7 n-l 


Modify : 


y n+l 


C n+ 1 + 


9 r l 

121 P n+1 c n+l 



Then the corrector may be iterated to convergence as desired. 
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E. P-C SETS CONSIDERED IN THE NUMERICAL EXPERIMENTS 



Having set forth the necessary foundation for the pre- 
dictor-corrector equations, a list of the P-C sets, which 
were analyzed and actually applied to the numerical solution 
of the initial value problem in (1-1) is now summarized. In 
each case the operator is evaluated in terms of y^ + ^, 

y ' i > ♦ • • 

7 n- 1 ’ 

P-C-I: Euler P-C set. 

The predictor is equation (2-7) and the corrector is 
(2-24), 



yjj+l = y n + hy^ ; T (x,h) = j h 2 y [2] (0 

\ 

■ v T (*- h > - - n> h V 31 «) 

P-C-II: Milne P-C set. 

Equation (2-10), in the case r = 3, q = 3 of the NBF 
extrapolation form, is used as the predictor, and the cor- 
rector is (2-16) where r = 1, q = 3, of the NBF interpolation 
form. This yields the equations 



n+1 

n+1 



>V3 * T 

y n -i * I 5 



T(x,h)« h 5 y [5] (Q 

T(x,h)=-|i y [5] (O 



P-C-III: Nystrom Predictor-Euler Corrector Set. 

The Predictor is (2-12) and the corrector is (2-24) which 
yields 

y l*i ' Vi + 2h K • T(x ’ h) ■ r y I31 (?) 
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Cl " y n + 7 ^A+l * y A> • T t x - 11 ) = - 17 C 131 CO 

P-C-IV: Hermite Predictor-Mi lne Corrector. 

The predictor used is (2-22) and the corrector is (2-16) 
yielding 

Cl " " 4y n + 5y n-l + h [ 4y A + ^A-l 1 

T(x.h) = ^ y' 4 > CC) 

0 

Cl ’ y n-l * 7 ly A-l * 4y A + y A+l ] 

TCx.h) = y 151 (?) 

P-C-V: Hamming Modified Predictor-Corrector Set. 

Equation (2-33) is the predictor and (2-34) is the cor- 
rector . 

Cl = y n-3 * r [2y A - y A-l + 2y A-2> 

T(x,h) - 11 h 5 y t51 CO 

Cl - I t 9y „- y „-2> * IT Iy A+l * 2y A - y A-l ] 

TCx.h) = C y [5, CO 

Adams Form (Adams-Bashforth Predictor-Adams -Moulton Corrector 
Set) . 

The A-B formulas are (in the case r = 0) equations of 
the NBF extrapolation form and the A-M formulas are (in the 
case r-0) equations of the NBF interpolation form. Three 
methods are considered: the second, third, and fourth order 

A-B predictor and A-M corrector. 
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P-C-VI: Second Order Adams P-C Set. 



The predictor and corrector are obtained in the case 
q=l yielding 

y n+l ' y n 4 T I3y n ' y A-l> 

T(x,h) = ^ h 3 y [31 (5) 

y n+ 1 ' y n + 7 Iy n + 1 + y A> 

3 

T (x , h) = ^ - y t3] (O 
12 

P-C-VII: Third Order Adams P-C Set. 

The case q = 2 in both the extrapolation and interpola- 
tion of the NBF form are used as the predictor and corrector 
equations respectively, resulting in the formulas 

y n+l = y n * TI I23y n * 16y n-l * 5y n-2 ] 

T(x,h) - i h 4 y l 41 «) 

y n+l = y n 1 * 8y n ' y n-l> 

Tfx.h) - ^ y [41 (O 

P-C-VIII: Fourth Order Adams P-C Set. 

The case q = 3 in both NBF forms yields the predictor- 
corrector equations 

y n+i - y n * rr [S5y A - 59y A-i * 37y A- 2 - 9y A-3> 
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T(x,h) 



251 

720 




CO 



y 



n+l 



= y 



n 24 



^A + i + 19 ^ - 5 ^-i + yA- 2 i 



T(x,h) = b 5y 5 CO 



From the P-C sets above an interesting fact is observed. 
Namely, P-C-I, P-C-III, and P-C-IV used different predictor 
equations but have the same corrector equation while P-C-II 
and P-C-V used the same predictor equation but have dif- 
ferent corrector equations. The significance of this 
observation will be clearly demonstrated later on. 
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V. NUMERICAL STABILITY OF PREDI CTOR- CORRECTOR METHODS 



In the numerical solution of an ODE, a sequence of 
approximations y n to the true solution y(x n ) is generated. 
Roughly speaking, the stability of a numerical method refers 
to the behavior of the difference or error, y - y(x ), as 
n becomes large. In order to begin the discussion, the 
various types of errors which are incurred in numerical 
integration of (2-4) must be considered. 

The errors incurred in a single integration step are 
of two types: 

1. The local truncation error or discretization 
error - the error introduced by the approximation of the 
differential equation by a difference equation. 

2. Errors due to the deviation of the numerical solu- 
tion from the exact theoretical solution of the difference 
equation. Included in this class are round-off errors, due 
to the inability of evaluating real numbers with infinite 
precision with the use of computer (i.e., computers usually 
operate on fixed word length) , and errors which are incurred 
if the difference equation is implicit and is not solved 
exactly at each step. 

If a multi-step method is used, an additional source of 
error results from the use of an auxiliary method (usually 
a single-step method, e.g., Runge-Kutta method), to develop 
the needed starting values for the multi-step method. 
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A. PROPAGATION OF ERROR IN A ONE-STEP METHOD 



The propagation of error in a one-step method can be 
analyzed by studying the particular representative one-step 
method of Euler (2-7). 



n+1 



= y + hf(x ,y ) 
1 n v n ’ 1 n' 



(2-36) 



This can be done by determining the relation of the error 
at step n+1 to the error at step n. To do this let 
denote the true solution of the initial value problem 
y ' (x) = f(x,y) with y(x Q ) = y . Then the total accumulated 
solution error e^ at step n is defined by 



The numerical values computed by Euler's algorithm 
satisfy the relation 



(2-37) 

(2-36) 



n+1 



= y 



n 



hf(x ,y ) - 
v n’ J n J 



R 



n+1 



(2-38) 



where R n+ ^ denotes the round-off errors resulting from 
evaluating (2-36). Similarly the true solution values 
satisfy the relation 



n+1 



= Y 



n 



hf <VV 



+ T 



n+1 



(2-39) 



where T n+ ^ denotes the local truncation error in (2-36). 
Subtracting (2-38) from (2-39) yields 

Y y = Y - y + h[f(x ,Y )- f(x ,y )]+E 
n+1 'n+1 n 'n 1 y n’ n' v n ,J n J n+1 

(2-40) 

where E n+ ^ = T^ + ^ + R n+ ^. Inserting (2-37) in (2-40) yields 
£ n.l ’ + h t f ( x n-V ' + E n + 1 C 2 ' 41 ’ 
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Applying the mean value theorem to the terms inside the 
bracket of (2-41), this relation between successive errors 
can be written as 

E n + 1 = Z n + h '(V>V f y Cx n>yJl * E nU 

where f denotes 9f/3y and y lies between y and Y . If 
y ' n J n n 

|fy(x,y)| _< C and |E n+ ^l £ E where C and E are both positive 
constants, the error expression above can be replaced by a 
related first-order difference equation (Henrici Ref. 5 
gives an excellent treatment of the solution of difference 
equations ) . 

e n +i = e n [ l + hC] + E (2-43) 

Under the stated assumptions it follows from (2-42) that 
|Z n+ ll £ |l n l [l*hC] .♦ E 

Now if |E I < e it follows that 
1 n 1 — n 

l E n + ll i e „ tl+hC] ♦ E 



Therefore, from (2-43) 

I E ■ , I < e A , 

1 n+1 1 — n+1 



(2-44) 



It follows then from (2-44) that 



\Z Q \ 1 e 0 I E i I < e i -*• 



.vrhere the symbol 
implies that. 



means 



That is, the propagated error is bounded by the solution of 

the related first-order difference equation. Now since 

Y -y„ = E„ = 0, the condition that |E I < e will be satis- 
o J o o ’ 1 o 1 — o 

fied by setting e Q = 0, which is the initial condition for 
the difference equation. 
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Defining G = 1+hC, the difference equation (2-43) can 



be written in the form 



e j , = Ge + E with initial condition e =0 
n+1 n o 

The solution of this first-order difference equation can 

be found by successive substitution as follows 



e l 

e 2 

e 3 



Ge o + 


E = 


E 


Ge-^ + 


E = 


(G+1)E 


Ge 2 + 


E = 


(G 2 +G+1)E 


Ge n- 1 


+ E 


= (G n_1 + G n " 2 + 



. . .+ G + 1)E 



The solution e n can be seen to be a geometric progression 
and as such can be written as 



e 

n 

It follows 
ithm (2-4) 



^G n - 1 
G - 1 



E 



\ / 

then that the propagated error 
is bounded by the expression 



in Euler's 



(2-45) 

algor- 



I E n+1 



(l+hC) n+1 -l 
hC 



We shall illustrate this estimate by a numerical example. 
Given the initial value problem 



y' = f(x,y) = -y with y(0) = 1 



we apply Euler's formula (2-36). It is required to express 
the error at x^ = 1 as a function of h assuming that the 

_ n 

round-off error |R n+ jJ £ 10 for single precision (double 
precision yields |R £ 10 ^) in the IBM 360/67 computer. 
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Computing 



C = max | f (x,y) | 

yields 

C = 1 

The truncation error associated with (2-36) is given by 

T , = - i h^y" (£ ) where x < £ < x 

n+1 2 1 '•^n-' n - n — n+1 

Knowledge of the second derivative y"(x) is needed to eval- 
uate this. Thus 

y"(x) = gl (y'fx)) - 4 (f(x,y)) 



= fx + fy 



dy 

dx 



Since f = 0 and f = -1, this reduces to 
x y ’ 

y"(x) = y. 



Using (2-45) , one obtains 



< (l+h) n -l 
n = h 



10 ^ + i h^ max | y | 
0<x< 1 



To evaluate max|y|, it is seen that the analytical solution 
of the initial value problem is 



Then for x = 1 
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* 















It is clear that max |y| is attained at x = 0 ; that is 
max |y| = 1. 

Substituting into the original expression yields 



jio -7 ♦ >d] 


(l+h) n - 1 


1 2 J 


h 



Rewriting this equation yields 



(2-46) 





j((l+h) n -l) 



(2-47) 



Thus the propagated error incurred in stepping forward from 
x = Otox= 1, as a function of h, is bounded by the above 
expression. 

The revised form of the propagated error bound shows the 
influence the various errors have in the numerical solution 
of the ODE. 

Rewriting (2-47) in the form 

e n = + r + 1 j((l+h) n -l) (2-48) 



where R ,, is still the round-off error 
n+ 1 

T^ + i is the truncation error reduced by 0(h) 



and ignoring the contents of the second parenthesis, the 
equation suggests that as h + 0, the truncation error 
decreases toward zero whereas the round-off error increases 
toward infinity. To minimize the error, it is clear that h 
must be chosen so that the truncation error and the round- 
off error have equal orders of magnitude. A numerical 
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experiment was conducted to verify this idea. Finally it 
can be seen that the accumulated error is not simply equal 
to the sum of the local truncation error and round-off error 
but must be computed as given by (2-48) which is the solu- 
tion of the first order difference equation corresponding to 
(2-36) . 

B. PROPAGATION OF ERROR IN THE MULTI-STEP METHOD 

As shown in equation (2-21) the most general representa- 
tion of the Predictor-corrector methods are of the form 



K K 

y., = X a.y .+h X B • y ' • 

'n+1 . ~ i 7 n-i . , l'n-i 

i = 0 i= - 1 



(2-49) 



However by confining the study to a form represented by 



K 

y = y + h X B-f • , f . = y' . 
'n+1 'n-r - , i n- i ’ n-i 'n-r 

i = - 1 



(2-50) 



simplicity is obtained without loss of generality. It can 
be verified that most of the formulas in Section III are in 
the form (2-50). For example, the case r = 0 and B_^ = 0 
the Adams -Bashforth formulas are obtained as 



K 

^n+1 ^n + ^ ®i^n-i 

while in the case r = 0 and B_^ f 0 the Adams -Moulton forms 
result . 

The propagation of error in the multi-step method can 
then be analyzed by studying equation (2-50). Following the 
same definition as in the previous section with regard to 
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, C, and E, 



Y , y , E , R , T , E 
n’ 7 n’ n’ n+1’ n+1’ n+1 

repetition can be avoided. 

The solution values generated by (2-50) satisfy the 
relation 

K 

>Vl ’ >Vr * 8 -l f °W>W + h J 0 Bi^n-i^n-d 



- R 



n+1 



(2-51) 



Similarly the solution values Y n satisfy the relation 

% 

Y n+1 - Y n-r + h ®-l £ (*n+l’ W + h j Q M <*n- i ’ Y n- i> 



+ T 



n+1 ‘ 



(2-52) 



Subtracting (2-51) from (2-52) yields the equation 



Y - y x 

n+1 7 n+l n-r 'n-r 



= _ t> - y^ _ ^ + he., Cf(x n+1 ,Y n+1 ) 



K 

- f (X ,y Al )) + h E 6 • (f (X . ,Y .) 
v n+l*'n+r ; • n i v n-i* n-i' 

i=0 



- f(x • ) + T - + R J ._ 
v n-i' n+1 n+1 



(2-53) 



By application of the mean value theorem to (2-53) and using 

the definition of E and E . . one obtains 

n n+1 



^n+1 ^n-r + ^-l^ Y n+l ^n+1^ £ ^ ^ x n+l ’^n+l^ 



(2-54) 



K 

+ h E 3 - (Y . -y .)fyCx • ,y -) + E ,, 
^_q n-i 'n-i' 7 v n-i’'n-i' n+1 

Using the definition of C,E, and E , (2-54) reduces to 



K 

E = E + hg ,E ±1 C + h E g. .C + E 
n+1 n-r -1 n+1 i n-i 
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or 



K 

E he n Z Al C = E + h I p.E .C + E 
n+1 -1 n+1 n-r ^_q i n-i 



(2-55) 



The related difference equation to (2-55) can be written as 

K 

e , (1 -h I 6 .|C) = e + hC E | 3 - | e . +E (2-56) 

n+i'- 1 -l 1 j n _ r i = o 1 n_1 

If I Z . I < e . (i = 0,1,... K) where K > r it follows that 
1 n-i 1 — n-i v ’’ J — 



K 



|s n+1 l( 1-hB.jCl) < |Z n . r | «■ hC S iBjl l^.jl * |E| 



Thus 



K 

l^ilCli-hB-iCl) 1 e-r + hC E I 6 

i=0 



e . + 
n - 1 



From (2-56) it follows that 

I ^ n +]_ I ( 1 1 1^ I ) .S e n+l (^”^l^-il) • 

If hC | B_ ^ | < 1, then 



(2-57) 



That is, if the magnitude of the propagated error is dom- 
inated by the solution of the related difference equation for 
K+l successive steps, namely |l^| £ e^ (i = 0,1,... K) then 
by induction, the same is true for all successive integer 
values of i. Therefore, a bound for the propagated error in 
the multi-step method can be determined by obtaining a 
solution of the related linear, inhomogeneous difference 
equation of (2.55). 
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As with ODE, (2-56) has homogeneous and particular 



solutions such that 



e = e rT + e 
n nH np 



The homogeneous solution is given by 



e nH = (l-hC|6_ 1 !)u K+1 -U K ' r - hC(|B 0 |u K 
+ | 6-lI u K_1 + ...+ |S K |) 



(2-58) 



The particular solution can be obtained by assuming that 
e n = -A, then substituting in (2-56) 



(1-hlB.JC) (-A)-(-A)-hC(-A) [|6 0 | + |B 1 | + .. . + |B k | 1 = E 



AhC E | B- | = E 
i= - 1 1 



hC E | B • | 
i=-l 1 

The general solution of the difference equation of (2-56) 
can be written as 



values and where u^ (i=l , 2 , . . . K+l) are the roots of the 
characteristic equation 



K 

A (hC | £ ,| + hC E |BJ] = E 
_i i=0 1 



K 




E 



(2-59) 



K 

hC E | B- 
i= - 1 1 
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In what follows the assumption is made that the roots are 
real and distinct (complex roots and multiple roots are 
discussed in detail in Ref. 4). 

It may be observed that the homogeneous part of the 
characteristic equation for the accumulated error (2-59) is 
exactly the same as in the characteristic equation for the 
original difference equation of (2-50). The significance 
of this point will be used in the analysis of numerical 
stability following a numerical example. 

To illustrate the previous analysis consider the case 
where r = 0 and K = 0 in (2-50) which yields 

y ., = y + £ [f X1 + f ) . (2-60) 

This is the modified Euler formula wherein y^ is given by 
the initial condition and the needed additional point has 
been supplied by a predictor formula. 

Assuming that f^(x,y) = -C and that E n+1 = E, the related 
difference equation yields 

(1 - | ( - C) ) u x - 1 - h ( -C) j = 0 
which has root 



U 1 = 



(i - ^) 

Cl + *y) 



(2-61) 



Using (2-59) the difference equation solution is then 



e n = d l 



h r -» n 
(i - ^) n 



(i + 



hC 
2 -J 



hC (i * i) 
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Analyzing the first term with as constant, it remains to 
show that if 




then the total error e n decreases with n and convergence is 
assured. Furthermore if this condition is met then 



and the solution of the difference equation decreases with 
increasing n. But the expression on the left of (2-62) is 
exactly the root of the characteristic equation of (2-60); 
thus it is clear that if u^ < 1, then the solution of the 
difference equation decreases with n. This idea can be 
extended to the general characteristic equation (2-50) whose 
solution may be written in the form 



It has been shown earlier that the total solution error 
I E . , I < e . , is bounded by the solution of the related 
difference equation. In turn it was observed that e n depends 
on the solution of the characteristic equation of the dif- 
ference equation. Therefore the problem of numerical sta- 
bility reduces to the analysis of the roots of the 
characteristic equation of the related difference equation 
of (2-63). 




(2-62) 



y = d,uV + d 0 u~ + 
'n 11 2 2 




(2-63) 
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C. STABILITY ANALYSIS OF A ONE-STEP METHOD 
Given an ODE 

y' = Ay with y(x Q ) = y Q (2-64) 

Assuming that A is constant, the analytic solution is given 
by 



y(x) = y Q e 



A < x n-V 



Using Euler formula (2-7), the related difference equation 
with y^ = Ay n is 

>Vi - ( 1 +Ah )y n - o 

whose characteristic root is 
u 2 = (1+Ah) 

Therefore the solution is 



uj = (1+Ah) n 



Since there is only one root, (2-63) reduces to 
7 n = Vj = d^ (1+Ah) n 



where d, = y . Therefore, 
1 ' o ’ 



y n = y 0 (l + Ah) n = y Q (1+Ah) 



We recall from calculus that 



e x - lim (l+h) x/h 
h 0 



(2-65) 
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and using the fact that A is constant (so that Ah •* 0 when 
h -*■ 0) , it can be seen that 



lim (1+Ah) 
h + 0 



A (x -x ) 
v n o J 

An 



It follows that (2-65) reduces 



A (x -x ) 
v n o J 
e 

to 



A (x 



n 



= y e 
1 o 



n 



■V 



The method is then (roughly speaking) stable, which means 
that a sequence of approximations y n will have the true 
solution as its limit if the corresponding values of h have 
zero as their limit. This analysis is quite simple since 
no extraneous roots represented by d^u^..., are i ntro ~ 

duced. The case where these roots are introduced will be 
discussed thoroughly in the next section. The only error 
introduced is represented by the particular solution of 
(2-59), whose behavior has been thoroughly studied in 
Section V-A as a function of h. 



D. STABILITY ANALYSIS OF A MULTI-STEP METHOD 
Given the same problem 

y' = Ay with y(x Q ) = y Q 

the corrector equation of Milne's P-C set (P-C-II) is used 
viz . 



n+1 



= y 



n-l 



+ h 

3 



y' • 

1 n - . 



4y * 
’ n 



y. 



n+l 



( 2 - 66 ) 



Since 




1 




1 
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yh = Ay n 



n+1 



= Ay 



n+1 



(2-66) has the related difference equation in the simplified 
form 



— 




"1 


1 




i— n 




~Ah 


! 




4Ah 


1, . Ah 


1 - 


3 _ 


1 


1' y n 


3 _ 


■ y n-l 1 + ~ 


_ 




J 






L J 



y n+i 

The characteristic equation yields 



u. 



with roots 



~i _ Ah " 


- Ui 






fi * 


3 _ 


1 


L 3 J 




3 



= 0 . 



(2-67) 



U 1 = 



2Ah ± / 9 + 5 Ah' 
3 - Ah 



Assuming that h is sufficiently small, the binomial theorem 
can be used to expand the square root and then simplifying, 
one obtains 

u^ = 1 + Ah + 0(h) 



U£ = -1 + -^ Ah + 0 (h) 



Since n = (x -x )/h, the solution of (2-67) can be written 



as 



y n = d^ (l+Ah+0 (h) ) 



(x -x )/h 

n 0 + d 2 (-l) n (l - i Ah + 0(h)) 






As h tends to zero, this solution tends to 



A(x -x ) -±A(x -x ) 

„ _ a ~ n o J . A ( ,-,n ^ 3 v n o' 

y = d 1 e + d 2 (-l) e 



( 2 - 68 ) 
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Since d^ is given by the initial condition and d£ is assumed 
to have been supplied by an auxiliary one step method, (2-68) 
can be completely determined. It is clear that there is an 
extraneous root introduced as a result of the use of a 2 n< ^ 
order difference equation to represent a first-order dif- 
ferential equation. The root u^ is then called spurious, 
parasitic, or extraneous and has no relation to the exact 
solution of the differential equation but, nevertheless, is 
unavoidable. This is clearly seen if it is assumed that 
d 2 = 0, in which case (2-68) reduces to 



y = y e 
7 n 7 o 



A(x -x ) 
v n o 7 



This is indeed the true solution of the ODE. But in general 
d 2 f 0, so the behavior of this extraneous root remains to 
be studied, since it affects the total solution. This could 
be done proving that, if A > 0, the parasitic solution 



~t( x _x ) 
n , ,.n 3 V n o 7 

u 9 = (-1) e 



tneds to zero exponentially as x n increases. Then the 
effect of the spurious solution becomes negligible as x n 
increases while the principal solution 



n 

1 



A ( x rT x o) 



tends toward the true solution values y (x ). However if 

n 

A < 0 the spurious solution increases exponentially in magni- 
tude and alternates in sign as the step-by-step calculation 
progresses, while the principal root u^ tends to one. 
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Moreover, at each stage of the calculation, rounding errors 
will introduce new spurious terms of the same type. Since 
these extraneous roots have no relation to the exact solu- 
tion, as they become large, they will dominate the results, 
thus invalidating the numerical solution. This situation is 
referred to as numerical instability. 

The same analysis can be generalized by going back to 
equation (2-63) 



y n d l U l + d 2 U 2 •” d K+l U K+l 



(2-69) 



n 



It is clear then that the principal root is d^u^ 



and 



, m j n 
d 2 2’ d 3 u 3 ’ 



n 



d K+l u K+l are extraneous roots that are intro- 



duced as a result of representing a first order differential 
equation by a (K+l) order difference equation. By analyz- 
ing the principal root u^ which has the form 

u^ = 1 + Ah + 0(h) 



or 

u^ = e^ + o (h) 

more precise requirements about numerical stability can be 
determined. Observe that u^ approximates e n ^ as was already 
shown . 

Comparing (2-69) to the exact solution of y* = Ay, viz, 
y(x n ) - e" Ah , 

it is clear that as long as | u^ | > |u^|, i = 2,3, ...K+l, 
when n' increases, the extraneous solution represented by u?, 
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* 



s- 





will become small when compared to the principal solution 
d^u^. In this situation the numerical solution is expected 
to be stable. However, if |u^| > | u^ | for any i = 2,3,...K+1, 
u 1 ? will become large with respect to u^ and the numerical 
solution component corresponding to u^ will dominate the 
solution. This situation is referred to as numerical 
instability . 

Intuitively then, the numerical solution will be valid 
if |ujJ > |u.|, i = 2,3, ...K+l. However, recall from the 
previous analysis of (2-59) that the characteristic equation 
(2-59) for the accumulated error e n+ ^, is exactly the same 
as the characteristic equation (2-63) for the original 
difference equation of (2-50). For a valid numerical solu- 
tion, it has been shown that the total error e must not 

n 

grow with n. From (2-59), it was observed that 

e nH = d l U l + d 2 U 2 + + d K+l U K+l 

and is equivalent to the condition that | u^ | <c 1, i , 2 , 3 , . . . K+l . 
It is then possible to define the stability of a method in 
the following way. 

(i) ABSOLUTELY STABLE if | u ± | <1, i = 1,2,... K+l 
(ii) RELATIVELY STABLE if | u ± | < u p i = 2,3, ...K+l 

Absolute stability does not imply relative stability. In 
other words, a numerical solution may have |u^| <_ 1, 
i = 1,2, ...K+l, but | u^ | < | u^| _< 1 , i = 2,3, ...K+l. 

Summarizing then, it can be seen that in the ODE 

y’ = Ay 
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absolute stability and relative stability are both important 
considerations if A < 0 since the exact solution is decreas- 
ing with x . If A > 0, the exact solution is growing with 
x n , so that relative stability is the important consideration. 
In other words, the numerical solution is valid as long as no 
component of u? dominates the component corresponding to the 
principal root. 



58 






* 

‘ 








VI. CONVERGENCE OF THE CORRECTOR IN THE P-C METHODS 



The term convergence has been mentioned earlier in the 

study of the propagation of the error as a function of h. An 

analysis of the convergence properties of correctors reveals 

that this convergence indeed depends on the step size h. A 

proper selection of h must therefore be made to ensure the 

convergence of the corrector. 

i t li 

Let y^ + i denote the j approximation of the solution 

value at x ... . Then the corrector formula can be written 
n+ 1 

in the form with f 0: 



. . K 

y** = y + h$ . f (x ,y J Al ) + h E £.f(x . ,y .) 

7 n+l •'n-r -1 v n+l’ ; n+l ; i=0 1 n-i ,7 n-i' 

(2-70) 

where i is the iteration index. Note that y 0 ., = y^ . , and 
J 'n+1 'n+1’ 

y-*ti is the corrector of iterating y-* , . which in turn is 
the corrector of (j-1) iteration. If the sequence of suc- 
cessive approximations y ^ , generated by iterating the 

* * 

corrector, converges to a limit, denoted by y then y n+ ^ 

satisfies the relation 



* * K 

^n+1 ^n-r + ^-l^ x n+l ’^n+1^ ^i^ X n-i ’^n- i^ 

(2-71) 

Subtracting (2-71) from (2-70) yields 



y 



^ ” y-n+-\ t^( X n+l ’^n+1-* ~ ^ ( X n + 1 > ^rTf 1 ^ ^ 



n+1 'n+1 



n+1 ’ 'n+1' 



which by the mean value theorem can be written 
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Cl ■ Ci ■ he -i [ Ci ■ Ci> Vvi'Vi' 

* i + 1 . , 

where y lies between y„., and y*: . If f < C, where 
/ n+l ' n+ 1 ' n+ 1 1 y 1 — 

* 

C is a positive constant, in the neighborhood of C^ n+ ^>y n+ ^)> 
then 



iCl - Cl> i h l e -ll C l y n+l ' Cll- 



It follows by induction that 



I C} - Cll 5 IhlB.alcjJ* 1 I y° tl - y* +1 | 



Therefore, the iteration of the corrector will converge to 
* 

the limit y n+ ^ provided that 

|he_ 1 C| < 1 (2-72) 

With this expression the conditions for convergence of the 
various corrector formulas considered can be written down 
immediately: 

„ Condition for 

Corrector _ - 1 Convergence 



Euler 


1/2 


hC 


< 


2 


Milne 


1/3 


hC 


< 


3 


Hamming 


3/8 


hC 


< 


8/3 


Third-order Adams 


5/12 


hC 


< 


12/5 


Fourth-order Adams 


3/8 


hC 


< 


8/3 



Both the Nystrom and second order Adams P-C sets used the 
Euler corrector while Hermite used the Milne corrector and 
thus are included in the above expressions. 
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An important point that must be brought out here is that 
convergence will only assure that the sequence of iterations 
y^i will converge to some definite value y , but not 
necessarily to the true solution y(x ). Stability is still 
the yardstick of a valid numerical solution though (2-72) 
provides a good rule of thumb to follow in the proper selec- 
tion of h. Indeed, later it will be shown experimentally 
that absolute stability can not be attained beyond the bound 
for h given by (2-72) since convergence is a necessary but 
not a sufficient condition for stability of a numerical 
solution (cf. Ralston [Ref. 6]). 
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VII. INFLUENCE OF THE PROPAGATED ERROR 



In order to test the idea presented in the previous sec- 
tion dealing with the propagation of error, the numerical 

example given in Section V-A was run on an IBM 360/67 compu- 

_ 7 

ter where the round-off err bound is 10 using single 
precision (SP) and 10 in double precision (DP). Recall 
that the ODE used was 

y' = -y with y(0) = 1 

and the final error expression using the Euler predictor 
formula at x = 1 is 

(l+h) n - lj (2-73) 

A series of h values from h = 10 ^ to h = 1.0 was used. The 
numerical values of e n were computed both in single precision 
and double precision. Table 1 shows the actual numerical 

_ 3 

data obtained. These data show that for large h to h * 10 , 

the SP and DP results are essentially identical, but as h 
decreases the round-off errors tend to increase while the 
truncation error tends to zero. The effect of the round-off 
error in DP remains negligible since, in double precision, 

14 digit accuracy is obtained. These results confirm the 
idea that the errors decrease with decreasing values of h 
to a certain point, beyond which the errors increase. A 
plot of the data obtained is shown in Figure 2. 
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TABLE 1 



INFLUENCE OF PROPAGATED ERROR (E) 
FOR y 1 =-y AT x=1.0 



h 


V * 

^ESP 


e edp 


1* 10' 6 


1.595-10" 1 


8. 763* 10' 7 


1-10' 5 • 


1 . 596 • 10 2 


8.593-10' 6 


5*10' 5 


3. 434 *10' 3 


4.295 -10' 5 


1*10" 4 


1. 780 *10' 3 


8.590 -10‘ 5 


5 * 10 ' 4 


9. Oil *10' 4 


2.147 *10" 4 


5 • 10 ~ 4 


7.722 *10' 4 


4.294 *10' 4 


1-10" 3 


1.029 *10" 3 


8.584 *10' 4 


5 • 10 " 3 


2.211 *10‘ 3 


2.143 -10‘ 3 


5 * 10 ~ 3 


4. 311 -10' 3 


4.278 *10‘ 3 


1-10" 2 


8.540-10' 3 


8. 524 *10' 3 


5 ' 10 " 2 


2. 106 *10' 2 


2.106 -10‘ 2 


5 * 10 " 2 


4.133-10' 2 


4.133 *10‘ 2 


1-10' 1 


7.968-10" 2 


7.968 -10 -2 


1*0 


5.0 -10" 1 


5.0 -10' 1 
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Error 



ESP - Error in Single Precision 
EDP - Error in Double Precision 




Figure 2. Influence of Propagated Error. 
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VIII. STABILITY BOUNDS OF P-C SETS 



Prior discussion about the stability of the predictor- 
corrector methods relates only to the limiting properties 
of the roots of the characteristic equation as the interval 
of integration approaches zero. In many applications, 
however, additional information about the actual size of h 
is needed to ensure stability of a method. The number of 
results discussed in published papers on this subject is 
extensive. But the interesting point observed is that 
although they all started with the analysis of the roots of 
the characteristic equations, they varied quite extensively 
in the mode of the predictor- corrector applications and in 
the number of P-C sets considered. To be more specific, 
some published papers are cited. The real negative stability 
expression used was hA where A = f in a single ODE. Chase 
[Ref. 7] analyzed the Milne P-C mode without iteration. The 
The real negative stability bounds found were - 0 . 8 < hA < - 0 . 3 , 
while on the other hand using the same P-C mode, but now 
iterated to convergence, resulted in numerical instabilities 
for all negative A. Hamming [Ref. 4] used three modes of 
his P-C set. First, iterated to convergence, the resulting 
stability bounds were -0.5 < hA < 0 ; second, with truncation 
error modification but without iteration the results were 
-0.85 < hA < 0; third, still with truncation error modifica- 
tion but with only two iterations, the results were 
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-0.9 < hA < 0 . Lapidus and Seinfeld [Ref. 2] using the PECE 
and PE(CE) S methods, where s is the iteration index, pub- 
lished the real negative stability bounds for all available 
predictor-corrector methods. Their algorithm is quite 
different from the algorithm used in this paper in that a 
final derivative evaluation is computed before terminating 
the computation with or without iteration. The symbol PECE 
or PE(CE) is used to denote their algorithms. The algorithm 
considered in this paper does not require a final derivative 
evaluation before termination with or without iteration, and 
hence can be denoted by P(EC) . 

To establish the real negative stability bounds for 
the P-C sets considered herein, the experimental procedure 
used by Chase to obtain the real negative stability bounds 
for hA will be followed. Chase used the test ODE 

y* = 100 - lOOy with y(0) = 0 

Using different values of h, he analyzed the behavior of the 

error until actual instability occurred. The P-C mode will 

be different in the sense that a standard application involv- 

2 

ing two iterations, P (EC) , will be used for all numerical 
methods considered. The choice of just two iterations was 
incluenced first by the published paper of Hull and Creemer 
[Ref. 8]. They conducted an experiment using the P(EC) 
mode applied to the Adams' methods only. After analyzing 
several numerical results, they concluded that the best 
method is S = 2, based on the cost of computation, average 
accuracy, and stability. Second, the choice of iterations 
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was influenced by the analysis of the corrector in terms of 
minimum truncation error. Consider Euler's P-C set (P-C-I), 
which has the form 



Predictor : 


y n+l 


II 

3 


+ hy ' 
1 n 




Corrector : 


y n + l 


= y n 


+ ll [y* 

2 iy n+l 


+ y'] , 
J n > 



as applied to the ODE y' = Xy. 

Ignoring the corrector, the single characteristic root 
is obtained from the predictor as follows: 

y - y - hXy = 0 
7 n+l 'n 'n 



u 1 - (1 + hX) = 0 
u^ = 1 + hX . 

Now apply the corrector once. Substituting the predictor 
into the corrector yields 

y„.i - d + >a ♦ y n 



The single characteristic root is 



u, = 1 + hX + 

1 2 



Now apply the corrector twice obtaining 



>Vi " y n + I dU * h* * y n * xy n i 



,2.2 . 3, 3 

y n+ i ■= [i ♦ h* * ♦ Lh_, y n 



The characteristic root is 
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1 + hX + 



U 1 = 



h 2 X 2 + h 3 X 3 



Continuing this process yields these results 



0 corrector: u, = 1 + hX 



1 corrector: u, = 1 + hX + 



hV 



2 correctors: u, = 



.2,2 .3,3 

1 + hX + 



3 correctors: u, = 



,2,2 ,3,3 ,4,4 

1 + hX + k A .... + ^ JL - + ^ -A 

2 4 8 



However, the true solution of y' = X^ is given by 

,, ,2,2 ,3,3 ,4,4 c 

e hX - 1 ♦ hX ♦ + M- ♦ ♦ °(h 5 ) 

2 o 24 

and thus the truncation error T(x,h) , given by 



u^ - e 



hX 



is, first with 0 corrector 

2 2 

T(x,h) = (1+hX) - (l+h\ + K^—+ 0 (h 3 ) ) 



T(x,h) = ^4^- - 0(h 3 ) 



By continuing this procedure with the 1,2,3 correctors and 
ignoring higher order terms, there results 



0 corrector: T(x,h) = 



-X 2 h 2 /2 



1 corrector: T(x,h) = 



-X 3 h 3 / 6 



2 correctors: T(x,h) = 

3 correctors: T(x,h) = 



X 3 h 3 /l 2 
X 3 h 3 /12 
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It can be seen that the use of the corrector more than twice 



seems to yield little advantage in this case. Thus the mode 
of P-C application to be used will be represented by the 
symbol P(EC) . In this mode of application the stability of 
the P-C set depends on both the predictor and corrector 
applications, though more so on the corrector equation. This 
assertion will be clearly illustrated in the case of the 
Milne, Hamming and Nystrom P-C sets. Thus, a priori know- 
ledge of the stability of a method can be obtained by 
analyzing the roots of the characteristic equation of the 
corrector formula. 

A. EXPERIMENTAL PROCEDURE 

Two general algorithms were developed to implement the 

predictor-corrector methods in Section IV-E. The first 

algorithm covers seven of the P-C sets, while the second 

t h 

algorithm covers the 8 P-C set, i.e., the Hamming modified 
predictor- corrector set. The flow chart for ALGORITHM 1 is 
on page 164 and that for ALGORITHM 2 is on page 167. How- 
ever, the general steps for each algorithm will be outlined 
here . 

Given the necessary starting values, the step size, the 
convergence term and the range of integration, the iterative 
predictor- corrector method for computing the numerical solu- 
tion y n+ i is set up as follows: 

ALGORITHM 1: 

Step 1: Compute the predicted value by the predictor 

formula. 
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Step 2 : 


Compute the derivative of the predicted value. 


Step 3: 


Compute the corrected value by the corrector 
formula . 


Step 4: 


Test for convergence. If convergence is attained 
go to step 7, otherwise proceed to the next step. 


Step 5 : 


Compute the derivative of the corrected value. 


Step 6: 


Compute the new corrected value by the corrector 
formula, using the value computed in step 3 as 
the new predicted value. 


Step 7 : 


The result is taken as the desired solution 
point. Advance the integration point by the step 
size. Return to step 1 to compute the next 
solution point. 



ALGORITHM 2: 



Step 1: 


Compute the predicted value by the predictor 
formula . 


Step 2: 


Modify the predicted value by adding the trun- 
cation error* of the predictor formula. 


Step 3: 


Compute the derivative of the modified predicted 
value . 


Step 4:. 


Compute the corrected value by the corrector 
formula . 


Step 5: 


Modify the corrected value by adding the trun- 
cation error* of the corrector formula. 


Step 6: 


Test for convergence. If convergence is attained 
go to step 10; otherwise proceed to the next 
s tep . 


Step 7: 


Compute the derivative of the modified corrected 
value . 



Step 8: Compute the new corrected value by the corrector 

formula using the value as computed in step 5 
as the new modified predicted value. 

* 

Truncation error as used in the formula is not the 
original truncation error but the modified truncation term 
expressed as a function of the difference between the cor- 
rected and predicted values. 
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Step 9: Modify the corrected value. 

Step 10: The result obtained is accepted as the desired 
solution point. Advance the integration point 
by the step size. Return to step 1 to compute 
the next solution point. 

2 

In both algorithms the mode of application is P(EC) as can 
be readily verified. In this scheme, if the corrector is 
used only once, then the number of function evaluations needed 
is also one, since it is assumed that the function value for 
the corrector has already been computed; if the corrector is 
used twice then two function evaluations are needed. Thus 
the number of function evaluations is determined by the 
number of iterations. 

The P-C algorithms were coded using the FORTRAN language. 
The fortran programs for P-C-I to P-C-VIII are listed on 
pages 170 to 200. In each program the meanings of symbols 
and parameters are explained through narrative comments. 

To determine whether the stability limits have been 
reached, certain criteria must be followed. Often, when the 
stability bound is approached through increasing step size h, 
the method begins to lose accuracy. An inaccurate, though 
stable, solution may often result in oscillations which 
appear at first to be due to actual numerical instability. 
Actual instability, associated with too large a value of h, 
usually results in an increasing oscillation in the error, 
or simply, the growth of the error in one direction. The 
particular error behavior depends on the sign of the para- 
sitic root causing the instability; a negative root causes 
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oscillations and a positive root causes a uniform error 
growth. 

B. PREDICTED STABILITY CHARACTERISTIC OF THE P-C SETS 

To gain a priori knowledge of the behavior of the sta- 
bility characteristics of the P-C sets, the roots of the 
characteristic equations of the corrector formulas will be 
studied. It must be noted that the behavior of the charac- 
teristic roots, though dominating the P-C method used as a 

set, will be influenced by the behavior of the predictor 

2 

formula. Since the chosen mode of application, P(EC) , is 
quite different from all the other modes of applications 
published in different papers, the stability limits obtained 
by several authors could not be used as the stability limits 
of the P-C set considered here. Thus the actual real nega- 
tive stability limits will be determined through numerical 
experiment. However, as noted previously, the analysis of 
the characteristic roots remains the same, and as such the 
published results are used and references are cited. This 
is done solely to conserve space and to avoid the tedious, 
repetitious process needed to solve the characteristic 
equations. The procedures have been amply presented and 
illustrated in detail in the section on numerical stability 
of predictor-corrector methods. 

a) Analysis of the characteristic roots of the corrector 
formulas . 

The P-C sets with common correctors will be grouped 
together except for the case of the second order Adams ' 



72 



method which has the same corrector as Euler, and which will 
be included in Adams' type. 

Euler (P-C-I) and Nystrom (P-C-III): 

The corrector formula is given by 

>Vl = y n + 7 + y n ] 

The related difference equation is in the form 

(1 - y n+1 - (1 + y n = 0, where C = f y (x,y) 

The single characteristic root is 

1 + llC 
2 

u i = ~ hC 
1 ' 2 

It can be seen that if C < 0, (i.e., the derivative 9f/3y is 
negative), then for absolute stability to occur, |u^| < 1 
which is satisfied if hC/2 < 1. Thus the numerical solution 
will decrease as does the exact solution If C > 0, as long 
as hC/2 < 1 then the numerical solution increases as does the 
exact solution. Figure 3 exhibits the behavior of the root 
versus hC as shown by Lapidus and Seinfeld [Ref. 2] . 

Milne (P-C-II) and Hermite (P-C-IV) : 

The corrector is of the form 

^n+1 ^n-1 + "3 ^n-1 + ^n + ^n+1^ 

The related difference equation is 



r, he - 




r 4hc“ 




^ he" 


r • -J 


>Vi - 


l— J 


^n - 
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-2.5 -2.0 -1.5 -1.0 -0.5 0 

hC 



Figure 3. Characteristic Root for Euler Corrector Formula. 
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The roots of the characteristic equation, as derived in 
section V-D, are of the form 



u^ = 1 + hC + 0(h) 
u 2 = -1 + i hC + 0(h) 



which can be written as 



u 



1 ~ 



e 



hC 



u 



2 “ 



-hC/3 

e 



If C > 0, u™ behaves like the exact solution and u^ dies out 
since, | u 2 | < 1. When, however, C < 0 , u^ decreases as 
does the exact solution, but u 2 increases. Thus the corrector 
is relatively stable but for C > 0, has no real negative 
stability bounds. Figure 4 shows the behavior of the roots 
as given by Chase [Ref. 7] . 



Hamming (P-C-V) : 

The corrector is given by 

>Vl " £ * X ^A+l + 2y n - 



The related difference equation is of the form 



1 3hC ' 




9 3hC 1 




I 3hC\ 


U - 1 ] 


y n+ l + 


8 4 J 


^n - 


1 8 J 



>Vl - 8 >V2 



= 0 



The characteristic equation is 



u 



. i + u 2 £ + 3hC 
8 1 U 8 4 



- u 



I 3hC 



= 0 



Chase [Ref. 7] analyzed the root behavior of this character- 
istic equation. Figure 5 shows the behavior of these 



1.6 




Figure 4. Characteristic Roots for Milne Corrector Formula. 
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1.6 



x - Positive Roots 
N - Negative Roots 
0 - Complex Roots 




Figure 5. Characteristic Roots for Hamming Corrector Formula. 
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characteristic roots versus hC. Analyzing the graph of 
Figure 5, it is evident that the corrector is stable for 
-2.4 <hC < 0. This corrector has also a wide range of 
relative stability as Chase has shown. In the graph the 
symbol used for positive roots is an x, that for negative 
roots is an N, and for the complex roots is an 0. 

ADAMS-MOULTON CORRECTORS (Second Order (P-C-VI, Third 
Order (P-C-VII), Fourth Order (P-C-VIII)): 

Since these formulas are all of the same type, they will 

be studied together and extended to the general case of the 

Adams -Moulton types. For the second order, the corrector is 

given by 

>Vl “ y n * 7 [ >W y; 1 - 

The related difference equation is of the form 

t 1 - I s ) >Vi - ii + y n = o 

with characteristic equation 

(1 - Jp] u - [1 + ££] = 0 

In the limit, as h -*■ 0, the root becomes 
u^ = 1 . 

For the third order, the corrector is 

Vl = ♦ 17 I5 ^ + l ♦ ■ ^-l 1 

The related difference equation is given by 
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0 



H - 12 hc ) >Vi - u * | hc ) y n + t§ y n -! ■ 

with characteristic equation 

[1 - jj hC] u 2 - [1 + | hC] u + = 0 



In the limiting case h -»■ 0 , the roots are 
u(u-l) = 0 
U 1 = 1 » u 2 = 0 • 



For the fourth order, the corrector is 

>vi = + It + 19 >A - S K- 1 * yA-2 1 - 

The related difference equation is 



[1 - 



3hC 



] y 



n+l 



- [1 



+ IT hC)y. 



n 



5hC 

24 



y 



n - 1 



-hC 

24 



y 



n- 2 



with characteristic equation 

f , 3 , 3 r i . 19 , 2 ^ 5hC hC _ n 

[1 " g hC] u - [1 24 hC]u 24 u ~ 24 ^ 



In the limiting case as h -> 0, the roots are 



u 2 (u-l) = 0 




Thus, it can be seen that, in each case, the dominant root 
is 1 while all the other roots (the parasitic ones) lie at 
the origin when h = 0. In the general q-step Adams -Moulton 
formula the equivalent characteristic equation is 



u q ' 1 (u-l) = 0 . 
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Thus, in every case, the A-M formulas have one root on the 
unit circle with all others at the origin when h = 0. 

As such they would seem to have excellent stability charac- 
teristics no matter what q-step formula is used. When h 
increases from zero, the parasitic roots move toward the 
unit circle. However, significantly large values of h can 
be reached before instability occurs. Crane and Klopfen- 
stein [Ref. 9] analyzed the behavior of the general Adams- 
Moulton formula and showed that it has stability bounds 
-1.3 < hC < 0 . It is interesting to note that, in general, 
as the order of the Adams corrector increases, accuracy 
increases but stability decreases. This is the main reason 
why higher order formulas are not considered (i.e., fifth, 
sixth, seventh, eighth). The graph of the characteristic 
roots of the general A-M formulas, as presented by Crane and 
Klopfstein, was not reproduced since it would give no 
additional information. 

In the foregoing analysis of the characteristic roots, 
the predicted stability limits served as a guide in the 
choice of the step size h to be used in the numerical 
experiment for obtaining the actual real negative stability 
bounds for each P-C set considered 

C. NUMERICAL RESULTS FOR REAL NEGATIVE STABILITY LIMITS 

In order to obtain actual real negative stability bounds 

? 

of the P-C sets considered in the P(EC) mode of application, 
the ODE 

y' = -y with y(0) = 1 
was chosen, where fy(x,y) = C = -1 < 0. 
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Each P-C set was run with increasing values of step 
size h, in small increments, to determine more precisely the 
actual stability limits. The criteria for actual numerical 
instability set forth previously, were followed in the 
analysis of data obtained. The prior knowledge of the roots 
of the characteristic equation obtained in the previous 
analysis narrowed the choice of starting values for the step 
size h for each method. Each method was run in single 
precision because of the large values of h. 

1. Euler CP-C-I) 

Since the predicted stability of the corrector is 
in the range -2 < hC < 0, the starting value chosen for h 
was 1.1, then was increased by 0.1 until instability occurred. 
Table 2 shows selected results for different values of h. 

At h = 1.1, it is absolutely stable since the error continues 
to decrease as x increases. The first sign of oscilla- 
tion was observed at h = 1.3 but it is still stable because 
the error is decreasing. Oscillation becomes evident at 
h = 1.6, but although the solution is inaccurate, it is 
not yet unstable, because there is no uniformly increasing 
trend of oscillation. Actual numerical instability occurred 
at h = 1.9 as is clearly shown by the uniformly increasing 
oscillation. Thus this method is stable within -1.9<hC< 0. 

As compared to the predicted stability region -2.0 < hC < 0 
obtained in [Ref. 2] where in the corrector is iterated to 
covergence, this was expected because of the use of only 
two iterations. 
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TABLE 2 



ABSOLUTE ERROR (E ) FOR EULER (P-C-I) 

ODE: y' = -y with y(0) = 1 

Single Precision 



h = 1.1 


h = 1.3 


E_ 




x E 


x E 


1.1 -2.239-10' 2 


1.3 -8.023-10' 2 


2.2 -5.765*10~ 2 


2.6 -1 . 019 * 10 _1 


3.3 -2.321*10' 2 


3.9 -2.748*10' 2 


4.4 - 1 . 6 10 * 10 ~ 2 


5.2 -2.957*10~ 2 


5.5 -6.078*10' 3 


6.5 -1.816‘10* 3 


6.6 -3.421*10' 3 


7.8 -8.301*10" 3 


7.7 -1.266-10' 3 


9.1 1.629*10' 3 


8.8 - 6 . 549 * 10* 4 


(STABLE) 


9.9 -2.406* 10," 4 




(STABLE) 




h = 1.6 


h = 1.9 


x Z E 


E„ 

x E 


1.6 -2.733*10 _1 


1.9 -6.696* 10 1 


3.2 -6.44* 10" 2 


3.8 5.308-10' 1 


4.8 -1 . 449 * 10 _1 


5.7 -1.546 


6.4 7 . 894 * 10 ~ 2 


7.6 3.020 


8.0 -1.442-10" 1 


9.5 -6.363 


9.6 1 . 6 34 * 1 0 ' 1 


(UNSTABLE) 


(STABLE BUT INACCURATE) 
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2. Milne (P-C-II) 

A priori knowledge indicates that this method has no 
real negative stability, so the starting value of h was 
chosen as 0.1. Table 3 shows the selected data obtained. 

Even at the starting value at h = 0.1, the buildup of error 
is consistent though small. Then at h = 0.4 oscillations 
become evident and the error tends to increase as x becomes 
large. At h = 0.7 the same increasing error tendency was 
observed. These observations, based on actual results, 
confirmed the predicted instability of tliis method for C <0. 
Thus it can be concluded that this method has no real 
negative stability bounds, as is evident from the actual 
data presented on Table 3. 

3. Nystrom (P-C-III) 

The selected starting value is close to the value 
used in Euler since this method has the same corrector. The 
program for Nystrom was run starting at h = 1.0 in increments 
of 0.1 until instability occurred. Table 4 shows the actual 
data obtained. The result for h = 1.1 to 1.3 show absolute 
stability. For h = 1.5, the error starts oscillating but not 
at a uniformly increasing rate. The solution corresponding 
to this data is inaccurate but stable. This inaccuracy 
continues up to h = 1.7. For h = 1.8 actual numerical 
instability has occurred as shown by the increasing oscilla- 
tion. Thus this method is stable for -1.8 < hC < 0. This 
range of stability is less than that for the Euler method 
though they have the same corrector equation. This confirmed 
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TABLE 3 



ABSOLUTE ERROR (Z) FOR MILNE (P-C-II) 

ODE: y'=-y with y(0)=l 

Single Precision 



h = 


0.1 


h = 0.4 


X 


I., 

M 


X 


M 


1.0 


no ' 7 


2.4 


5 . 1 7 * 10 ~ 5 


2.0 


5 * 10 " 7 


3.6 


-5.86’10' 5 


3.0 


5* 10' 7 


4.8 


1.115‘10" 4 


4.0 


5 * 10 " 7 


6.0 


-1.555-10' 4 


5.0 


6 * 10 ~ 7 


7.2 


2 . 393* 10' 4 


6.0 


1-10" 6 


8.4 


-3.58-10' 4 


7.0 


1.6'10" 6 


9.6 


5 . 394 * 10 " 4 


8.0 


2. 1- 10" 6 




(UNSTABLE) 


9.0 


3.0*10' 6 






10.0 


4 . 2 * 10 6 






(UNSTABLE) 






h = 


0.7 






X 








2.8 


2 . 5 84 ’ 10 ~ 4 






4.2 


7 . 106 ’ 10~ 4 






5.6 


8. 312* 10" 4 






7.0 


1. 186-10' 3 






8.4 


1.918-10' 3 






9.8 3.213-10' 3 

(UNSTABLE) 
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TABLE 4 



ABSOLUTE ERROR (X) FOR NYSTROM (P-C-III) 

ODE: y' = -y with y(0)=l 

Single Precision 



h = 1.1 


h = 1.3 


X 


e n 


X 




2.2 


3.465-10* 2 


2.6 


5.570-10* 2 


3.3 


' 4. 369*10' 2 


3.9 


5. 976* 10" 2 


4.4 


3.043-10' 2 


5.2 


5.013-10' 2 


5.5 


1. 724*10" 2 


6.5 


1 . 446 * 10' 2 


6.6 


7.181-10' 3 


7.8 


1 . 914 ‘ 10" 2 


7.7 


3. 062 * 10~ 3 


9.1 


-1. 183* 10‘ 2 


8.8 


4. 488* 10* 4 




(STABLE) 


9.9 


3 . 363 * 10 " 4 








(STABLE) 






h = 1.5 


h = 1.8 


X 


Z N 


X 


S N 


3.0 


7.005-10' 2 


3.6 


-1.865* 10' 2 


4.5 


6. 538* 10' 2 


5.4 


2.057*10' 1 


6.0 


8. 581 ' 10' 2 


7.2 


-3. 134* 10' 1 


7.5 


- 2 . 784 ’ 10 ' 2 


9.0 


1.132 


9.0 


1 . 128 * 10 1 


10.8 


-3.2066 


(STABLE BUT INACCURATE) 




(UNSTABLE) 
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the idea that the stability of the predictor- corrector 
depends on both the predictor and corrector formula. 

4. Hermite (P-C-IV) 

This method has the same corrector as the Milne 
method, hence the choice of starting value is the same, 
h = 0.1, which is then incremented by 0.1 until instability 
occurred. Table 5 shows the selected actual data obtained. 
From these results it is evident that this method is unstable 
since for all values of h the error growth is steadily 
increasing. This is expected since it was predicted that 
the corrector has a dominant role in the stability of this 
method as a P-C set. Thus this method, which uses a Milne 
corrector, followed the instability behavior of the Milne 
corrector, within two iterations. Therefore, it can be 
seen that the Hermite P-C set has no real negative stability 
bounds . 

5 . Hamming (P-C-V) 

The chosen starting value was h = 0.5 with increment 
0.1. The predicted stability of the corrector is -2.4 < hC < 0. 
Table 6 shows, that for h = 0.5 to h = 0.6, the method is 
absolutely stable. For h = 0.7 to h = 0.8 the method is 
stable although the numerical solution becomes inaccurate. 

For h = 0.9 the error grows in one direction thus resulting 
in actual instability. Thus the method is stable within 
the range of -0.9 < hC < 0. Compared to the predicted sta- 
bility of the Hamming corrector, given by -2.4 < hC < 0, the 
range of stability is greatly reduced. This is due to the 
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TABLE 5 



ABSOLUTE ERROR (E) FOR HERMITE (P-C-IV) 

ODE: y ' = -y with y(0)=l 

Single Precision 



h - 


0.2 




h = 


= 0.5 


X 


Z HE 


X 




y 

^HE 


1.0 


-1*10“ 6 


1.0 




2.820-10' 4 


2.0 


6 * 10 " 6 


2.0 




3. 966* 10' 4 


3.0 


- 5 . 1 * 10 " 6 


3.0 




6. 146- 10' 4 


4.0 


9.1-10' 6 


4.0 




9. 347* 10" 4 


O 

LO 


-1.18-10' 5 


5.0 




1 . 401 ‘ 10" 3 


6.0 


1 . 71 ' 10 5 


6.0 




2.088*10' 3 


7.0 


- 2 . 4 * 10 " 5 


7.0 




3 . 106 ‘ 10” 3 


8.0 


3 . 39 * 10 " 5 


8.0 




4.617*10" 3 


9.0 


-4.78-10' 5 


9.0 




6.863*10' 3 


10.0 


6 . 74 * 10 " 5 


10.0 




1.019*10' 2 


(UNSTABLE) 




(UNSTABLE) 


h = 


* 0.8 








X 


V 

^HE 








1.6 


2.388-10' 3 








3.2 


6 . 331 ‘ 10 " 3 








4.8 


1 . 521 * 10 2 








6.4 


3. 298* 10" 2 








8.0 


7.033*10 -2 








9.6 


1 . 496 * 10 _1 








(UNSTABLE) 
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TABLE 6 



ABSOLUTE ERROR (E) FOR HAMMING (P-C-V) 

ODE: y' = -y with y(0)=l 

Single Precision 



h = 0.5 


h = 0.6 


x Z HA 


x Z HA 


2.0 -2 . 056 ' 10~ 4 

3.0 -1 . 418 * 10" 4 

4.0 -8.06'10" 5 

5.0 -4. 15’10 -5 

6.0 -2.02*10" 5 

7.0 -0 . 6 * 10 " 6 

8.0 - 4 . 6 ’ 10 -6 

9.0 -2.2’10 -6 

10.0 -l.l’lO' 6 

(STABLE) 


2.4 -3. 741’10~ 4 

3.6 - 2 . 897 * 10 _4 

4.8 -1.689*10 4 

6.0 -9.37'10~ 5 

7.2 - 5 . 62 * 10” 5 

8.4 -3.78*10" 5 

9.6 - 2 . 81 ’ 10" 5 
(STABLE) 


h = 0.7 


h = 0.8 


h = 0.9 


x E HA 


x Z HA 


x Z HA 


2.8 - 5 . 709 ’ 10" 4 
4.2 -5.717’10 4 
5.6 - 3 . 944 ' 10" 4 
7.0 - 3 . 04 7 ’ 10 ” 4 
8.4 -2.874-10' 4 

9.8 - 3. 001 ' 10~ 4 
(STABLE) 


3.2 -7.459*10 -4 

4.8 -1.096’10' 3 

6.4 -9. 78 ’ 10 ” 4 

8.0 - 1 . 1 21 * 10 " 3 

9.6 -1 . 534* 10" 3 

(STABLE BUT 
INACCURATE) 


3.6 -8.327-10 -4 

5.4 -1 . 982 ’ 10" 3 

7.2 -2.335’ 10" 3 

9.0 -3.801*10~ 3 

(UNSTABLE) 




effect of the Milne -predictor . But recalling the discussion 

of the published results, this actual real negative stability 

bound is exactly the same as that obtained by Hamming using 

two iterations with truncation error modification. This 

2 

coincidence is not without basis since the P(EC) mode has 
exactly the same truncation error modification as that of 
the Hamming method, making the two experimental procedures 
identical . 

6 . Second Order Adams (P-C-VI) 

The starting value chosen was h = 0.1 with an 
increment of 0.1. Table 7 shows the actual d^ta obtained. 

For h = 0.1 up to h = 0.3 the method is stable. But for 
h = 0.4 the method shows actual instability as evidenced by 
the uniformly increasing error growth. Thus this method 
is stable within the range -0.4 < hC < 0. Compared to the 
predicted stability of the corrector, which was the same as 
Euler CP-C-I] , with range of stability -2.0 < hC < 0, it is 
evident that the stability is greatly reduced. Again this 
could be traced to the high instability of the predictor 
formula used. This is probably the reason why almost all 
published material on the Adams-Moulton method always starts 
with the third order and higher predictor. Hull and Creemer 
[Ref. 8] showed that this method has the largest error among 
the Adams -type methods. 

7 . Third Order Adams (P-C-VII) 

The choice for h starts at h = 0.8 and is then 
incremented by 0.1. Table 8 shows the selected actual data 
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* 






TABLE 7 



ABSOLUTE ERROR (E) FOR SECOND ORDER ADAMS (P-C-VI) 

ODE: y ' = -y with y(0)=l 

Single Precision 





h = 


0. 


1 


O 

II 

r£, 


. 2 


X 






E A2 


X 


E A2 


1.0 






3. 017- 10' 4 


1.0 


1.130- 10' 3 


2.0 






2. 862- 10' 4 


2.0 


1.347- 10' 3 


3.0 






2 . 223- 10' 4 


3.0 


1.277- 10' 3 


4.0 






1. 772* 10' 4 


4.0 


1. 196- 10' 3 


S.O 






1 .527- 10' 4 


5.0 


1.146-10' 3 


6.0 






1.408- 10' 4 


6.0 


1.120-10~ 3 


7.0 






1. 354- 10" 4 


7.0 


1. 108- 10' 3 








-4 




, -3 


8.0 






1.330-10 


8.0 


1.103-10 


9.0 






1.320-10" 4 


9.0 


1.100-10' 3 


10.0 






1. 315- 10' 4 


10.0 


1.099- 10’ 3 




(STABLE) 


(STABLE) 




h = 


0. 


3 


h = 0 


.4 


X 






S A2 


X 


e A2 


1.2 






2.607-10’ 3 


1.2 


4.164-10' 3 


2.4 






3.S89- 10' 3 


2.4 


7. 290- 10' 3 


3.6 






3. 814- 10' 3 


3.6 


8 . 663 • 10 " 3 


4.8 






3. 861- 10“ 3 


4.8 


9. 207- 10' 3 


6.0 






3. 865- 10' 3 


6.0 


9 . 411 • 10 " 3 


7.2 






3. 869- 10 _3 


7.2 


9.484- 10‘ 3 


8.4 






3. 868- 10' 3 


8.4 


9. 510- 10' 3 


9.6 






3. 868- 10' 3 


9.6 


9.519- 10' 3 


(STABLE 


BUT INACCURATE) 


(UNSTABLE) 
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TABLE 8 



ABSOLUTE ERROR (£) FOR THIRD ORDER ADAMS (P-C-VII) 

ODE: y’=-y with y(0)=l 

Single Precision 



h = 0.8 


h = 0.9 


X 


Z A3 


X 


Z A3 


3.2 


•1.402-10" 2 


3.6 


1 . 600 • 10~ 2 


4.8 


7. 751- 10" 3 


5.4 


8.272-10" 3 


6.4 


1. 968* 10' 3 


7.2 


1.832* 10" 3 


8.0 


4.614-10' 4 


9.0 


7.436* 10' 4 


9.6 


4.608* 10' 4 




(STABLE) 




(STABLE) 






h = 1.2 


h = 1.3 


X 


E A3 


X 


E A3 


3.6 


7.06-10' 5 


3.9 


-4. 736* 10' 3 


4.8 


8. 773*10' 3 


5.2 


1.426* 10' 3 


6.0 


5. 019* 10' 3 


6.5 


1. 569* 10' 3 


7.2 


5. 216- 10' 3 


7.8 


1 . 869 * 10" 3 


8.4 


-2.976-10' 3 




(UNSTABLE) 


9.6 


3 . 050 * 10" 3 






(STABLE BUT INACCURATE) 
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obtained. For h = 0.8 to h = 0.9 absolute stability occurs. 
For h = 1.0 to h = 1.2 oscillation occurred and the solution 
is inaccurate though stable. For h = 1.3 the error growth 
is increasing showing actual instability. Thus the method 
is stable within the range -1.3 < hC < 0, which is in 
conformity with the predicted stability of the general Adams- 
method, -1.3 < hC < 0. This also justifies why most Adams- 
Moulton methods start with the use of the third order predic- 
tor. 

8 . Fourth Order Adams (P-C-VIII) 

Starting value was h = 0.5 with 0.1 increment. 

Table 9 shows the actual data obtained. From h = 0.5 to 
h = 0.7 the method is absolutely stable. For h = 0.9 oscil- 
lation continues as with h = 0.8. This behavior was observed 
up to h = 1.0 but the numerical solution is still stable 
though inaccurate. For h = 1.1 actual instability is evident 
from the tendency to increasing oscillation of the error. 

Thus this method is stable for -1.1 < hC < 0. 

Summarizing then, the real negative stability bounds of 
the different P-C sets considered are outlined below with 
C = fy(x,y) and h = step size. 



P-C Set 


Real Negative 
Stability Limit 


P-C- 1 


-1.9 < hC < 0 


P-C-II 


Unstable 


P-C-III 


-1.. 8 < hC < 0 


P-C-IV 


Unstable 


P-C-V 


-0.9 < hC < 0 


P-C-VI 


-0.4 < hC < 0 


P-C-VII 


- 1.3 < hC < 0 


P-C-VIII 


-1.1 < hC < 0 
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TABLE 9 



ABSOLUTE ERROR (E) FOR FOURTH ORDER ADAMS (P-C-VIII) 

ODE: y' = -y with y(0)=l 

Single Precision 



h = 0.5 


h = 0. 7 


X 


E A4 


X 


e A4 


2.0 


5 . 6 7 • 1 0 " 5 


2.8 


4.094-10" 4 


3.0 


2.125-10' 4 


4.2 


8. 157- 10‘ 4 


4.0 


1.488-10" 4 


5.6 


3. 522- 10' 4 


5.0 


8. 11 • 10‘ 5 


7.0 


1. 239* 10" 4 


6.0 


3.95-10" 5 


8.4 


3. 52- 10' 5 


7.0 


1.8-10 ‘ 5 


9.8 


7. 7- 10" 6 


8.0 


7.9- 10‘ 6 




(STABLE) 


9.0 


3. 4- 10' 6 






10.0 


1.4-10' 6 








(STABLE) 






h = 0.9 


h = 1.1 


X 


e A4 


X 


e A5 


3.6 


1. 819- 10' 3 


3. 3 


-3. 885’ 10 3 


5.4 


2 . Ill* 10 3 


4.4 


5. 789’ 10' 3 


7.2 


4 . 93 • 10 " 5 


5.5 


1 . 214’ 10' 3 


9.0 


-5.125-10' 4 


6.6 


3. 818’10' 3 


(STABLE BUT INACCURATE) 


7.7 


6. 321’ 10‘ 3 






8.8 


-4.534* 10' 3 






9.9 


7. 846* 10' 3 








(UNSTABLE) 
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The significant points noted from these results are: 

i) Within the limits of their stability, the Hamming 
and Fourth Order Adams methods produced the best accuracy. 

ii) It is quite interesting to see that although the 
Euler and Nystrom methods have the widest range of real 
negative stability, their results are not quite as accurae 
as the results of the Hamming and the Fourth Order Adams, 
even the Third Order Adams methods. 

D. NUMERICAL RESULTS FOR RELATIVE STABILITY 

If C > 0, the solution itself and generally also the 
error are increasing exponentially. In this situation 
relative stability is the important consideration. A P-C 
method is relatively stable if the rate of change of the 
error with respect to a finite range of integration points 
is less than the rate of change of the true solution with 
respect to the same finite range of integration. By this 
definition it could be seen that it is extremely difficult 
to establish a fixed bound for relative stability. In 
order to have a working knowledge of the relative stability 
of the P-C sets considered the ODE 

y 1 = y with y(0) = 1 

where C = fy(x,y) > 0 was used. Each P-C set was run with 
the same series of h values from h = 0.1 to h = 2.0, with 
an increment of 0.1 The range of integration for each 
step size h was x = 0 to x = 10. Table 10 shows the true 
solution values from x = 1 to x = 10.0. Table 10-A to 
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TABLE 10 



TRUE SOLUTION VALUES FOR 

ODE: y'=y with y(0)=l 

Single Precision 



X 


^exact 


1.0 


2.718 


2.0 


7.389 


3.0 


20.085 


4.0 


54.597 


5.0 


148.409 


6.0 


403.417 


7.0 


1096.595 


8.0 


2980.838 


9.0 


8102.710 


10.0 


22025.320 
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Table 10-H show the selected actual data obtained for the 
different P-C sets. 

Table 10-A shows that the Euler method is relatively 
stable for the particular range of integration. Comparing 
the error growth with the solution growth of Table 10, it 
could be seen that the rate of change of error is less than 
the rate of change of the solution. From h= 0.1 toh= 0.5 
the method is quite accurate. However from h = 1.0 to 
h * 2.0 it becomes quite inaccurate though still relatively 
stable. Table 10-B shows that the Milne method is not only 
relatively stable but accurate from h = 0.1 to h = 2.0. 

From Table 10-C the Nystrom method exhibits relative stabi- 
lity from h = 0.1 to h = 1.0 but becomes unstable at h = 2.0. 
This could be determined by simply comparing the magnitude 
of the error at x = 10.0 and h = 2.0 with the magnitude of 
the true solution at the same point, in which case it is 
obvious that the error is larger than the true solution. 

This is another way of looking for relative instability 
since it is clear that if the rate of change of the error 
is greater than the rate of change of the solution with 
respect to n (the number of solution points) then eventually 
the error will be greater than the true solution. Table 10-D 
Shows that the Hermite method is relatively stable and also 
accurate. From Table 10-E it can be seen that Hamming's 
method is also stable and accurate. Table 10-F shows that 
the second order Adams' method is stable from h = 0.1 to 
h = 1.0, but the accuracy is quite diminished at h = 1.0. 

At h = 2.0 it is unstable. From Table 10-G, the third 
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TABLE 10 -A 



ABSOLUTE ERROR (E) FOR EULER (P-C-I) 

ODE: y ' =y with y(0)=l 

Single Precision 



v 


e e 




h = 0.1 


h = 0.5 


h = 1.0 


h = 2.0 


1.0 


-2.231- 10' 3 


-4 . 819 * 10~ 2 


-1.567* 10 1 


- 


2.0 


-1.219-10' 2 


-2.570-10 -2 


-7.515-10 -1 


-1.610 


3.0 


- 4 . 980 ’ 10~ 2 


-1.046 


-2.959 


- 


4.0 


-1 . 806 * 10 _1 


-3.808 


-10.638 


-16.401 


5.0 


-6.142-10' 1 


-13.011 


-36.261 


- 


6.0 


-2.004 


-42.721 


-119.355 


-155.571 


7.0 


-6.360 


-136.450 


-383.286 


- 


8.0 


-19.774 


-427.070 


-1208.460 


-1420.042 


9.0 


-60.496 


-1316.117 


-3756.500 


- 


10.0 


-182.828 


-4006.574 


-11546.150 


-12622.530 
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TABLE 10 -B 



ABSOLUTE ERROR (E) FOR MILNE (P-C-II) 

ODE: y'=y with y(0)=l 

Single Precision 





Z M 


A 


h = 0.1 


h = 0.5 


h = 1.0 


h = 2.0 


1.0 


4. 8- 10' 6 


- 


- 


- 


2.0 


-2.38-10" 5 


1.409* 10' 3 


- 


- 


3.0 


-1.678-10" 4 


-1.831* 10' 4 


- 


- 


4.0 


-7.477-10' 4 


-1. 309* 10' 2 


3.070* 10' 1 


- 


5.0 


-2.777-10" 3 


-7. 122* 10' 2 


8. 342* 10' 1 


- 


6.0 


-9. 765’ 10' 3 


-2 .917- 10" 1 


1.794 


- 


7.0 


-3. 222- 10‘ 2 


-1.059 


4.026 


- 


8.0 


-1.040*10‘ 1 


' -3.608 


8.289 


362.921 


9.0 


-3.281'10 _1 


-11.781 


15.589 


- 


10.0 


-1.003 


-37.402 


23.242 


3472.707 
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TABLE 10 -C 



ABSOLUTE ERROR (E) FOR NYSTROM (P-C-III) 

ODE: y'=y with y(0)=l 

Single Precision 







X 


h = 0.1 


h = 0.5 


h = 1.0 


h = 2.0 


1.0 


-2.034*10' 3 


-2. 756- 10' 2 


- 


- 


2.0 


-1. 173- 10" 2 


-2.231-10' 1 


-5. 244* 10 -1 


- 


3.0 


-4.875-10 -2 


-1.016 


-2.758 


- 


4.0 


-1. 783* 10" 1 


-3.901 


-11.304 


-16.401 


5.0 


-6 . 096 * 10" 1 


-13.757 


-41.718 


- 


6.0 


-1.996 


-46.134 


-145.106 


-215.571 


7.0 


-6.350 


-149.630 


-485.911 


- 


8.0 


-19.778 


-473.892 


-1584.737 


-2384.042 


9.0 


-60.589 


-1474.339 


-5069.113 


- 


10. 0 


-183.289 


-4523.753 


-15975.800 


-24472.530 
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TABLE 10 -D 



ABSOLUTE ERROR (£) FOR HERMITE (P-C-IV) 

ODE: y'=y with y(0)=l 

Single Precision 





Z HE 


X 


h = 0.1 


h = 0.5 


h = 1.0 


h = 2.0 


1.0 


3. 8* 10" 6 


-3.948-10" 4 


- 


- 


2.0 


-2.77-10" 5 


-2.104-10" 3 


-1.217-10" 2 


- 


3.0 


-1.831-10" 4 


-9.292-10" 3 


-2.285-10' 2 


- 


4.0 


-7.629-10" 4 


-3.546-10" 2 


-8. 180- 10" 2 


7.092-10" 1 


5.0 


-2.853-10" 3 


-1. 244- 10" 1 


-2.519-10" 1 


- 


6.0 


-9. 765‘ 10' 3 


-4. 147- 10" 1 


-7. 766- 10 _1 


21.334 


7.0 


-3. 247-10" 2 


-1.335 ' 


-2.355 


- 


8.0 


-1.044-10 -1 


-4.194 


-7.-67 


275.2517 


9.0 


-3.281-10" 1 


-12.941 


-21.023 


- 


10.0 


-1.003 


-39.367 


-62.078 


2860.671 
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TABLE 10-E 



ABSOLUTE ERROR (E) FOR HAMMING (P-C-V) 

ODE: y'=y with y(0)=l 

Single Precision 





Z HA 


A 


h = 0.1 


h = 0.5 


h = 1.0 


h = 2.0 


1.0 


2.19-10" 5 


- 


- 


- 


2.0 


6 . 77 ' 10" 5 


2.980-10' 3 


- 


- 


3.0 


1 . 8 31 • 1 0 ” 5 


2 . 79 2 * 10 ' 3 


- 


- 


4.0 


6.714-10' 4 


-7.202-10’ 3 


4 . 332 ' 10* 1 


- 


5.0 


1.831-10" 3 


-6.004- 10" 2 


8.498- 10' 1 


- 


6.0 


5. 371* 10 -3 


-2. 777- 10" 1 


1.043 


- 


7.0 


1 . 879 • 10” 2 


-1.039 


-5. 24- 10' 1 


- 


8.0 


4.687' 10' 2 


-3.639 


-10.532 


380.836 


9.0 


1.406* 10" 1 


-12.097 


-53.558 


- 


10.0 


4.296-10' 1 


-38.882 


-213.566 


3040.449 
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TABLE 10 -F 



ABSOLUTE ERROR (£) FOR SECOND ORDER ADAMS (P-C-VI) 

ODE: y '=y with y(0)=l 

Single Precision 





Z E 


X 


h = 0.1 


h = 0.5 


h = 1.0 


h = 2.0 


1.0 


-2.100-10' 3 


-2. 735*10‘ 2 


- 


- 


2.0 


-1.242*10‘ 2 


-2. 358* 10' 1 


-5.041-10' 1 


- 


3.0 


-5.241-10" 2 


-1.127 


-2.792 


- 


4.0 


-1.934-10' 1 


-4.453 


-11.936 


-14.401 


5.0 


-6.649* 10" 1 


-16.013 


-45.345 


- 


6.0 


-2.186 


-54.444 


-161.08 


-201.571 


7.0 


-6.974 


-178.434 


-548.354 


- 


8.0 


-21.766 


-569.834 


-1812.745 


-2324.042 


9.0 


-66.800 


-1785.156 


-5866.652 


- 


10.0 


-202.378 


-5510.332 


-18684.180 


-24506.530 



102 



TABLE 10 -G 



ABSOLUTE ERROR (E) FOR THIRD ORDER ADAMS (P-C-VII) 

ODE: y'=y with y(0)=l 

Single Precision 



Y 


Z A3 


A 


h = 0 . 1 


h = 0.5 


h - 1.0 


h = 2.0 


1.0 


2.155-10' 4 


- 


- 


- 


2.0 


1.318-10" 3 


2 .076-10' 2 


- 


- 


3.0 


5.569-10' 3 


1. 454- 10" 1 


2.131-10' 2 


- 


4.0 


2.062-10" 2 


6 . 352 ’ 10 ~ 1 


5.538-10' 1 


- 


5.0 


7.075-10' 2 


2.376 


2.827 


- 


6.0 


2 . 324 * lO" 1 


8.216 


11.147 


14,938 


7.0 


7.412-10' 1 


27.090 


39.645 


- 


8.0 


2.307 


86.509 


132.954 


217.016 


9.0 


7.089 


269.992 


429.273 


- 


10.0 


21.468 


828.183 


1349.750 


2409.128 
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TABLE 10 -H 



ABSOLUTE ERROR (Z) FOR FOURTH ORDER ADAMS (P-C-VIII) 

ODE: y'=y with yCO^l 

Single Precision 





Z A4 


A 


h = 0.1 


h = 0.5 


h = 1.0 


h = 2.0 


1.0 


8 . 6 ‘ 10 " 6 


- 


- 


- 


2.0 


-2.77-10" 5 


5 . 5 3 * 10 ' 5 


- 


- 


3.0 


-1.984'10‘ 4 


-1. 995* 10‘ 2 


- 


- 


4.0 


-8.087-10" 4 


-1. 089‘ 10 _1 


9.992- 10' 2 


- 


5.0 


-3.234-10' 3 


-4 . 450 * 10 _1 


-7. 852- 10" 1 


- 


6.0 


-1 .074* 10‘ 2 


-1.614 


-5.094 


- 


7.0 


-3. 540-10' 2 


-5.490 


-21.956 


- 


8.0 


-1. 176- 10" 1 


-17.926 


-81.859 


262.678 


9.0 


-3.476-10' 1 


-56.878 


-283.230 


- 


10.0 


-1.085 


-176.812 


-936.164 


1748.47 
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order Adams' method showed accuracy and stability from 
h = 0.1 to h = 2.0. Finally, Table 10-H shows that the 
Fourth order Adams' method is quite stable and accurate 
from h = 0.1 to h = 2.0. 

The most interesting points observed from these actual 
results are: 

i) Though the Milne and Hermite methods have no real 
negative stability limits, they are quite accurate for 

C > 0. In fact the Hermite method produced the least start- 
ing error, and the Milne's method provided the second least 
starting error, and both maintained accuracy up to h = 2.0. 

ii) In contrast, the Euler and Nystrom methods both 
have the widest range of real negative stability limits but 
showed inaccurate results starting at h = 1.0. In fact the 
Nystrom method is unstable at h = 2.0. 

iii) The third and fourth order Adams' and Hamming methods 
showed wide range of relative stability and provided accurate 
results. The second order Adams' method again showed inac- 
curate results at h = 1.0 and relative instability at h = 2.0. 

iv) From h = 0.1 to h = 1.0, the methods in terms of 

accuracy, rank as follows: Milne ranks first, Hermite second, 

then Hamming, Fourth Order Adams, and Third Order Adams in 
that order. 
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IX. NUMERICAL EXPERIMENTS 



In order to test the performance of the different P-C 
sets considered, a collection of test ODEs with many unusual 
and interesting features (i.e., singularities, discontinui- 
ties, infinite derivatives, oscillating derivatives, etc.) 
were selected. Several of such ODEs were presented by Hull 
and Creemer [Ref. 8] and Lapidus and Seinfeld [Ref. 2] . 

Table 11 lists the test ODEs considered. To simplify nota- 
tion, the ODE number will be used to specify the differential 
equation. For example ODE I is equivalent to the ODE 
y» = _y + I0sin3x with y(0) = -3, whose analytic solution is 
y(x) = sin3x - 3cos3x. 

From the previous section on stability bounds of P-C 
sets it was observed that the different methods showed good 
accuracy up to h = 0.5 on both experimental ODEs, y' = -y 
and y' = y, considered. It was further observed from pre- 
vious numerical results that all the P-C sets exhibited 
good stability behavior within the range of h up to h = 0.5 
except for the Milne and Hermite methods, in the case of 
real negative stability. But these shortcomings of the Milne 
and Hermite methods were compensated by the fact that they 
produced good accuracy and wide range of relative stability. 
This analysis is needed in the sense that the test ODEs con- 
sidered exhibit both cases of fy(x,y) < 0 and fy(x,y) > 0. 

By using values of h = a ^ , 2 ^ , 2 ^ , 2 ^ , 2 ^ , 2 ^ , 2 ^ 
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where the largest value of h = 0.5, each method will be 
expected to perform quite well. Thus a comparative analysis 
of their accuracy and computing time would be meaningful. 

In dealing with accuracy, the effects of the propagation 
of errors would be strongly felt, as has already been shown 
in the previous analysis of error propagation. While 
computing time takes into its fold the number and complexity 
of function evaluations. 

All the eight P-C sets were run on each test ODE using 
the series of h values as previously mentioned. The Fortran 
programs used were the same as those used for stability 
analysis, but the precision used was Double Precision (14 
digit precision for IBM 360/67 computer) to minimize the 
effect of rounding error as h decreases. Each test ODE will 
be studied individually, with the objective of providing 
useful and important recommendations on which P-C set 
performs best for the particular class of problems. 

A. ODE I 

Table 12 shows the exact solution from x = 1.0 to 
x = 10.0. From Table 12-A the results show that the Euler 
method provides the best accuracy with 7.65 secs in computing 
time with the Nystrom's ranking second in accuracy though 
with less computing time, 7.46 secs. The Milne and the 
Hermite methods are not considered since they both showed 
instability in their numerical solution values. This could 
easily be seen by looking at Table 12-B which shows the 
behavior of the errors for these methods. The errors for 
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TEST ODEs 
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TABLE 12 



TRUE SOLUTION VALUES FOR ODE I 
Double Precision 



X 


^exact 


1.0 


3.111097 


2.0 


-3.159926 


3.0 


3.145509 


4.0 


-3.068134 


5.0 


2.929351 


6.0 


-2.731937 


7.0 


2.479843 


CO 

o 


-2.178115 


9.0 


1.832792 


10.0 


-1.450785 
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TABLE 12 -A 



ODE I AT x = 10.0 

ABSOLUTE ERROR FOR P-C-I, P-C-II, P-C-III, and P-C-III 

Double Precision 



1/h 


Z E 




e n 


r 

^HE 


128 


-9. 83479- 10' 5 


4.51483- 10' 9 


-9.83522-10" 5 


-8. 53288* 10' 3 


64 


-3.93360-10' 4 


8.21731* 10' 8 


- 3 . 93415 ’ 10 " 4 


- 3. 40492 ’ 10" 2 


32 


-1. 57286- 10' 3 


-1. 17807* 10" 7 


- 1. 57372* 10' 3 


-1.35592 * 10 _1 


16 


- 6 . 2 8126 * 10 ~ 3 


-2. 51544- 10‘ 5 


-6 . 92664 • 10' 3 


-5. 38364-10' 1 


8 


-2.49436-10' 2 


-1. 35324- 10' 3 


-2. 52237- 10' 2 


-2.134035 


4 


-9 . 64382 * 10" 2 


-3.35324-10' 2 


- 1. 00786 * 10" 1 


-6.039341 


2 


-3. 34570-10' 1 


1.527752 


-3. 38212* 10 _1 


-37.089294 


CPTS* 


7.65 


5.59 


7.46 


8.37 



CPTS means computing time in seconds. 
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TABLE 12 -B 



ABSOLUTE ERROR FOR P-C-I, P-C-II, P-C-III, and P-C-IV 

Double Precision 
ODE I ; h = 0 . 5 



X 


e e 






y 

^HE 


1.0 


7.171-10" 1 


- 


4. 159* 10 _1 


-1.584 


2.0 


-4 . 640 • 10 _1 


2.081-10" 1 


-6. 649* 10' 1 


-1.645 


3.0 


5. 756- 10" 1 


9.447*10" 2 


5.817-10' 1 


-2.464 


4.0 


-5 . 392 * 10" 1 


2 . 509 * 10 _1 


-6.091-10" 1 


-3.386 


5.0 


5.418-10" 1 


2.076-10" 1 


5. 806* 10 _1 


-5.218 


6.0 


-5.154-10' 1 


4 . 2 66 * 10 * 1 


-5.560-10" 1 


-7.537 


7.0 


4 . 854 ’ 10~ 1 


4.737-10 -1 


5 . 148* 10 -1 


-11.400 


8.0 


-4.431* 10' 1 


7.910’10 _1 


-4.652* 10' 1 


-11.747 


9.0 


3.929-10' 1 


1.011 


4. 056' 10' 1 


-25.066 


10.0 


-3.345 


1.527 


- 3 . 382 ’ 10" 1 


-37.089 
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the Milne and Hermite methods increased in one direction 
eventually growing to magnitudes greater than the true 
solution at x = 10.0 (i.e., for Milne's error, 1.527 > 1.450 
and for Hermite's error -37.089 >1.450), thus these methods 
are not suited for solving ODE I. The Euler and Nystrom 
methods both started at accuracy approximately equal to 
fifth decimal places and ended up at h = 0.5 with accuracy 
up to the first decimal place. But though their accuracy 
is quite restrained as h increases, nevertheless the behavior 
of the error tends to decrease as the range of integration 
is increased. From this group then the Euler method is 
the best choice for ODE I. 

In Table 12-C the results for the other four methods are 
listed. All these methods showed stability for all values 
of h used, but Hamming's produced the best accuracy and 
the best computing time. Hamming's method started at 10 ^ 
accuracy at 1/h = 128 and ended with 10 at h = 0.5. The 
Fourth Order Adams' comes next in both accuracy and computing 
time, then the Third Order Adams', and the Second Order 
Adams' came in that order. 

Comparing these last four methods with the Euler method, 
only Adams' second order method is inferior. Therefore, for 
solving ODEs that belong to the class of ODE I the Hamming 
method is highly recommended in terms of accuracy and least 
cost in computer time. The Fourth Order Adams' is the next 
choice in the absence of the Hamming method. 
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TABLE 12 -C 



ODE I at x = 10.0 

ABSOLUTE ERROR (E) FOR P-C-V, P-C-VI, P-C-VII, and P-C-VIII 

Double Precision 



1/h 


^ Ha 


Z A2 


E A3 


Z A4 


128 


3.0381-10' 10 


-9.86551-10' 5 


-9. 60790' 10" 7 


1 . 860 70 * 10 " 8 


64 


1 . 12626 • 10 ~ 8 


-3.95790-10' 4 


-7. 81167* 10 -6 


3.18596-10* 7 


32 


4 . 68494 * 10 " 7 


-1.59200-10" 3 


-6. 45309- 10' 5 


4.04497-10' 6 


16 


1 . 01544 * 10 5 


-6 . 42959 ’ 10~ 3 


-5.49889-10' 4 


5.91417-10" 5 


8 


3.69762* 10' 4 


- 2 . 60 335 * 10 ~ 2 


-4.97320- 10' 3 


7. 48403' 10" 4 


4 


1.29180-10" 2 


-1. 03170' 10' 1 


-5.00046'10' 2 


4.56055'10~ 3 


2 


6.39906- 10' 2 


-3.586413 


-4. 74203- 10' 1 


-1. 90291' 10 _1 


CPTS 


5.76 


7.98 


8.55 


5.99 
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B. ODE II 



The true solution values are given in Table 13 from 
x = 1.0 to x = 10.0. From the data presented in Table 13-A 
it is clear that the most important consideration is the 
step size h. For h = 0.5 the four methods identically 
yield errors greater than the true solution at the same 
point. In solving problems of this type then h must be 
chosen to be smaller than 0.25 for the numerical solution 
to be valid using the Euler, Milne, and Nystrom methods. 

The Hermite method obviously is not worth considering, 
because of its inaccurate results for all values of h. 

For h <0.25 the Milne method gives the best numerical 
solution and the least computing time, followed by the Euler 
then the Nystrom methods. 

Table 13-B shows that to have a vlid numerical solution 
for any of the other four methods h must be less than 0.25. 
For values of h < 0.25 , the Hamming method yields the 
greatest accuracy, followed by the Fourth, Second, and 
Third Order Adams'. The Milne method is a little better 
than the Second Order Adams'. Thus for accuracy the Hamming 
method ranks first, then the Fourth Order Adams, Milne, 

Euler, Second Order Adams, Third Order Adams, and lastly 
the Nystrom methods. 

For problems in the class of ODE it is recommended that 
h must be chosen less than 0.25 or smaller if valid numerical 
solution and accuracy are desired. Then use Hamming as the 
numerical method to solve the ODE. Again the Fourth Order 
Adams method is the second best choice. The choice of h can 
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TABLE 13 



TRUE SOLUTION VALUES FOR ODE II 
Double Precision 



X 


^exact 


1.0 


- 1.381773 


' 2.0 


- 4 . 931505 * 10 _1 


3.0 


8 . 488724 * 10 _1 


4.0 


1.410446 


5.0 


6 . 752620 - 10" 1 


6.0 


- 6 . 807547 • 10 -1 


7.0 


- 1.410888 


8.0 


- 8.438582 * 10" 1 


9.0 


4 . 990117 ’ 10" 1 


10.0 


1.383092 
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TABLE 13-A 



ODE II at x = 10.0 

ABSOLUTE ERROR FOR P-C-I, P-C-II, P-C-III, and P-C-IV 

Double Precision 



1/h 


Z E 






Z HE 


128 


2.596512-10' 3 


2 . 50786 ' 10" 9 


8.444203*10' 4 


-4.446340-10' 1 


64 


2.518638-10' 3 


1.634916* 10' 7 


7.950519-10" 3 


-1.764591 


32 


2 . 726949 * 10~ 3 


9.921730-10' 6 


5.455366-10' 2 


-6.947361 


16 


4.095807-10' 2 


5.69316710' 4 


4 . 234564 • 10 _1 


-26.911798 


8 


6.213302 -10' 1 


2 . 178410 * 10" 3 


3.214692 


-100.834152 


4 


9.160642 


1. 362692* 10 _1 


23.407199 


-353.009679 


2 


131.604044 


5.482274 


157.762746 


-1081.991963 


CPIS 


6.38 


4.69 


5.34 


8.26 
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TABLE 13-B 



ODE at x = 10.0 

ABSOLUTE ERROR FOR P-C-V, P-C-VI, P-C-VII, and P-C-VIII 

Double Precision 



1/h 


e ha 


E A2 


Z A3 


Z A4 


128 


-1. 20446* 10 9 


-3. 760418-10' 4 


4. 448046* 10' 4 


-2. 334506* 10" 8 


64 


-7. 48373* 10' 9 


-3. 192068* 10' 3 


3. 610822* 10' 3 


-6 . 282905 * 10 7 


32 


5. 531148* 10' 7 


-2 . 331998 * 10" 2 


2 .961989* 10' 2 


-1. 347495* 10' 5 


16 


7. 353901* 10' 5 


-1.532622-10' 1 


2.455737-10' 1 


-6 . 597818 * 10' 5 


8 


3.686296* 10' 3 


-7.673571-10' 1 


2.016912 


-1.600424-10' 2 


4 


1. 657253* 10' 1 


-6. 706531* 10" 1 


15 . 053093 


-2. 268373* 10' 1 


2 


6.292332 


44.871662 


80.187459 


1.401846 


OPTS 


5.26 


7.79 


8.36 


4.96 
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be formally stated as |hC| < 0.25, since the value of h to 
be used is directly dependent upon C = f^,(x,y). 

C. ODE III 

Table 14 shows the exact solution values for ODE III from 
x = 1.0 to x = 10.0. From the data of Table 14-A it is clear 
that all four methods performed quite well to a good degree 
of accuracy for all balues of h. But for best accuracy and 
least computing time the Milne Method stands out followed by 
the Euler, Hermite, and Nystrom methods in that order. 

From Table 14-B, again it is observed that any of the last 
four methods yield a valid numerical solution, to a certain 
degree of accuracy. For an excellent accuracy, however, the 
Hamming method is the best choice. Its range of accuracy is 
from 10 up to 10 ^ for values of h ranging from 1/h = 128 
to 1/h = 2. The Milne method compared to these last four 
methods ranks second only to the Hamming, though its com- 
puting time is a little smaller than Hamming's. Thus for 
excellent accuracy, the order of choice is an follows: 

Hamming, Milne, Fourth Order Adams, Third Order Adams, Her- 
mite, Euler, Second Order Adams, and lastly the Nystrom 
method in solving ODEs belonging to the class of ODE III. 

D. ODE IV 

True solution values are listed in Table 15 from x - 1.0 
to x = 10.0. Actual data obtained in Table 15-A showed that 
all four P-C sets are good numerical methods for ODE IV. But 
if a choice is to be made the Milne's accuracy is far greater 
than the other three methods with computing time only 0.16 
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TABLE 14 



TRUE SOLUTION VALUES FOR ODE III 
Double Precision 



X 


^exact 


1.0 


8. 414709* 10 _1 


2.0 


9 . 092974 * 10 _1 


3.0 


1. 411200* 10 -1 


4.0 


-7. 568024* 10 _1 


5.0 


-9.589242* 10 _1 


o 

vo 


-2.794154* 10 _1 


* 

o 


6. 569865* 10 _1 


o 

• 

00 


9 . 893582* 10 _1 


9.0 


4. 121184* 10 _1 


10.0 


-5. 440211* 10 -1 
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TABLE 14 -A 



ABSOLUTE ERROR FOR P-C-I, P-C-II, P-C-III, and P-C-IV 

ODE III AT x=10 . 0 
Double Precision 



1/h 


e e 


S M 


e n 


e he 


128 


-1. 705737*10' 6 


-1. 099* 10" 11 


-1. 785639-10" 6 


1. 341400* 10" 7 


64 


-6.845247*10"^ 


9 . 286 * 10" 11 


-7.814645-10" 6 


1.288789-10" 6 


32 


-2 . 703738-10" 5 


1. 26041* lO 9 


-3.531438-10" 5 


1.031052 *10' 5 


16 


-1.089843-10" 4 


1. 58770* 10" 8 


-1. 736454* 10' 4 


8.247845-10' 5 


8 


-4.365210-10' 4 


1. 596440* 10' 6 


-9. 528314* 10' 4 


6 . 594418 * 10" 4 


4 


-1. 755574-10' 3 


4. 317133* 10' 5 


-5.856889-10" 3 


5. 260049 *10" 3 


2 


-7.174787-10' 3 


1 . 071599 * 10" 3 


-3.910558* 10' 2 


4 .154229* 10" 2 


CPTS 


7.81 


5.59 


7.12 


7.92 
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TABLE 14 -B 



ODE AT x = 10.0 

ABSOLUTE ERROR FOR P-C-V, P-C-VI, P-C-VII, and P-C-VIII 

Double Precision 



1/h 


V 

nA 


l AZ 


Z A3 


E A4 


128 


2.04* 10' 12 


5. 023475- 10' 6 


-1.900550- 10" 8 


3.683* 10" 11 


64 


6. 396- 10" 11 


-7.798238- 10' 6 


-6.720151- 10' 7 


6 . 5237 ’ 10~ 10 


32 


2.01370* 10" 9 


-3. 512798- 10" 5 


-5. 373282* 10' 6 


1. 246830- 10' 8 


16 


6. 236927* 10' 8 


-1. 721024- 10' 4 


-4.287630- 10' 5 


2. 649671- 10' 7 


8 


1. 970409* 10' 6 


-9.408691- 10‘ 4 


-3. 391988- 10‘ 4 


6.946177-10' 6 


, 4 


6. 368642- 10' 5 


- 5 . 767408 * 10~ 3 


-2.588546- 10' 3 


1. 805833- 10' 4 


2 


1. 569286*10' 3 


-3. 848313- 10‘ 2 


-1.686135- 10* 2 


4. 360227- 10' 3 


CPTS 


6.23 


10.0 


10.0 


8.23 
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TABLE 15 



TRUE SOLUTION VALUES FOR ODE IV 
Double Precision 



X 


^exact 


1.0 


3.894003 


2.0 


3.639183 


3.0 


3.306565 


4.0 


2.943035 


5.0 


2.578543 


6.0 


2.231301 


7.0 


1.911513 


8.0 


1.624023 


9.0 


1.370189 


10.0 


1.149189 
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TABLE 15 -A 



ODE IV at x = 10.0 

ABSOLUTE ERROR FOR P-C-I, P-C-II, P-C-III, and P-C-IV 

Double Precision 



1/h 


£ e 


M 


Z N 


V 

l he 


128 


-7.344900-10' 7 


-3. 4* 10' 3 


-4. 598339- 10" 7 


4 . 575814 • 10" 6 


64 


-1.839350-10" 6 


-6.07-10' 12 


-1.840047-10' 6 


1. 826764- 10" 5 


32 


-7. 335362-10" 6 


-1. 0603* 10" 10 


-7. 339871*10' 6 


7. 278726- 10" 5 


16 


-2.926381-10" 5 


-1.97060-10' 9 


-2. 872267' 10" 5 


2.889080-10' 4 


8 


- 1 . 176812 * 10" 4 


-4.063982-10' 8 


-1.119506-10" 4 


1. 128099* 10" 3 


4 


-4. 709013- 10' 4 


-6. 024899- 10' 7 


-4. 260573- 10" 4 


4 . 418111 * 10" 3 


2 


-1. 885932* 10' 3 


-5.791009*10"^ 


-1.541727-10' 3 


1.668370-10' 2 


OPTS 


4.39 


3.83 


3.67 


6.74 
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sec longer than the next best choice which is Nystrom. After 

Nystrom, the Euler method is preferred, then the Hermite. 

Table 15-B data again demonstrates that any of the last 

four P-C sets can be used for solving ODE IV. Milne method 

compared with these methods competes with the Hamming in 

accuracy and has an edge in computer time (3.83 secs < 4.45 

secs). But this fraction of seconds difference is overcome 

by the demonstrated better accuracy of the Hamming in every 

- 7 - 1 

step of the solution from h = 2 to h = 2 . Thus, if 

preference is to be made, Hamming is the most likely choice 
as the best method for solving ODE IV class of problems. 

The Milne, Fourth Order Adams, Third Order Adams, Nystrom, 
Second Order Adams, Euler, and Hermite are the order of 
choices following the Hamming method. 

E. ODE V 

Table 16 shows the true solution values for the range of 

integration considered. From the data of Tables 16-A and 

16-B, it is obvious that any method can be used to solve 

ODE V and obtained accuracy up to minimum of 10 ^ and 

maximum of 10 Thus in comparing these methods, accuracy 

criteria must not be the greatest concern. By studying 

closely the behavior of the errors as h increases, it was 

noted that the Milne, Hamming, and Fourth Order Adams 

-7 -5 

exhibited decreasing errors from h = 2 to h = 2 for 
Milne, and from h = 2 ^ to h = 2 ^ for both Hamming and 
Fourth Order Adams. Analyzing further the Hamming and 
Fourth Order Adams methods it could be seen that the Fourth 
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TABLE 15 -B 



ODE IV at x = 10.0 

ABSOLUTE ERROR FOR P-C-V, P-C-VI, P-C-VII, and P-C-VIII 

Double Precision 



1/h 


Z HA 


E A2 


Z A3 


E A4 


128 


- 1 3 

3.0*10 i:> 


-4.585327* 10" 7 


1.28499*10' 9 


-1 . 61 ’ lO 12 


64 


- 13 

2.4*10 1-5 


- 1. 829584* 10 _6 


1. 025906-10' 8 


-2.690* 10" 11 


32 


- 1 2 

2.37*10 


-7.282200*10' 6 


8.174763-10" 8 


-4 .3747* 10' 10 


16 


7.849-10' 11 


-2.884380*10' 5 


6 . 488701 ’ 10" 7 


-7. 20867* 10' 9 


8 


2. 78944*10' 9 


-1 . 131593* 10' 4 


5. 111696* 10' 6 


-1.225493*10' 7 


4 


1.186192-10" 7 


-4. 359162* 10 4 


3. 970077 * 10" 5 


1 . 778276* 10' 6 


2 


2 .202890 * 10' 6 


-1.625329*10^ 


3.006535* 10' 4 


2.20 1592 * 10 ' 5 


CPTS 


4.45 


6.46 


7.02 


4.02 
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TABLE 16 



TRUE SOLUTION VALUES FOR ODE V 
Double Precision 



X 


v 

1 exact 


1.0 


1 . 020204 


2.0 


1.040832 


3.0 


1.061913 


4.0 


1.083472 


5.0 


1.105541 


6.0 


1.128152 


7.0 


1.151338 


8.0 


1.175139 


9.0 


1.199593 


10.0 


1.224744 
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TABLE 16 -A 



ODE V at x = 10.0 

ABSOLUTE ERROR FOR P-C-I, P-C-II, P-C-III, and P-C-IV 

Double Precision 



1/h 


Z E 


M 


e n 


Z HE 


128 


2.26254' 10‘ 9 


7.0' 10" 14 


-1.94002'10" 9 


-4.8'10" 13 


64 


9.04072-10" 9 


4. 0' 10' 14 


-7.747767'10" 8 


-4.45' 10" 12 


32 


3. 609166' 10" 8 


0* 


-3. 088771' 10' 8 


-3.583' 10" 11 


16 


3.436058’10' 8 


-2.2*10 i:> 


-1. 227293-10" 7 


-2.8689-10" 10 


8 


4 . 960094 ’ 10 " 7 


-3. 28' 10' 12 


-4.844143-10' 7 


-2.29641'10" 9 


4 


1.978142' 10" 6 


-4 . 237 ’ 10" 11 


-1.886736' 10' 6 


-1.839847'10" 8 


2 


7.866820'10' 6 


-4.0408' 10‘ 10 


- 7 . 644850 ’ 10" 6 


- 1 . 478586 ' 10" 7 


CPTS 


5.63 


5.58 


5.16 


5.35 



* 

0 means error less than 10 
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TABLE 16- B 



ODE V at x = 10.0 

ABSOLUTE ERROR FOR P-C-V, P-C-VI, P-C-VII, and P-C-VIII 

Double Precision 



1/h 


Z HA 


Z A2 


Z A3 


J A4 


128 


4.0-10' 13 


-9.29021*10" 9 


-6. 0* 10" 14 


1.6-10' 13 


64 


1.9-10' 13 


-1. 179555* 10' 8 


-2.46* 10' 12 


1.1-10' 13 


32 


1 . 2 * 1 0 " 1 3 


-3. 331471* 10" 8 


-2.018- 10" 11 


7 . 0 * 10 " 14 


16 


-14 

5.0*10 


-1.249903-10" 7 


-1.6114*10" 10 


-3.0*10' 14 


8 


- 1 3 

-1.4*10 1-5 


-4. 926236* 10' 7 


-1.27940*10‘ 9 


- 1 . 14 * 10" 12 


4 


-4. 83* 10" 12 


-1.950214* 10' 6 


-1.007531-10" 8 


- 1 . 802 * 10* 11 


2 


-1.1744-10' 10 


-7.649028 *10' 6 


-7. 798880 * 10~ 8 


-2. 7389-10" 10 


OPTS 


6.31 


10.0 


10.0 


7.78 
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- 4 

Order Adams had its optimum value of-h at h < 2 while that 

- 3 

of Hamming achieved its optimum h value at h < 2 . Thus 

the Hamming method has a higher value of h for its range of 
excellent accuracy. This property gives the Hamming method 
a slight edge over that of Milne, and the Fourth Order Adams. 
However if both accuracy and computing time are considered 
the Milne method is likely to be the choice. Therefore for 
the class of ODE V, the recommendation for the best method 
to be used will be most likely dependent upon the particular 
interest: whether accuracy and wider range of h values are 

the primary concern, or whether accuracy and least computing 
time is the criterion. For the former criteria the Hamming 
method serves the best purpose while for the latter the 
Milne method offers the best solution. However, as was 
noted earlier, if no criterion is involved but the interest 
is just to solve the problem the easiest way, the self- 
starting Euler method is recommended. For realistic pur- 
poses however the following order of choice is highly 
recommended: Hamming, Milne, Fourth Order Adams, Third 

Order Adams, Hermite, Nystrom, Euler, and lastly Second 
Order Adams . 

F. ODE VI 

True solution values are listed in Table 17 from x = 1.0 
to x = 10.0. Results from Table 17-A clearly showed that h 
must be chosen to be quite small to have a valid numerical 
solution. For the series of h values chosen none of the 
methods exhibited valid numerical solution from 1/h = 32 up 
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TABLE 17 



TRUE SOLUTION VALUES FOR ODE VI 
Double Precision 



X 


^exact 


1.0 


1.732050 


2.0 


2.236067 


3.0 


2.645751 


4.0 


3.0 


5.0 


3.316624 


6.0 


3.605551 


7.0 


3.872983 


O 

oo 


4.123105 


9.0 


4.358898 


10.0 


4.582575 
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TABLE 17- A 



ODE VI at x = 10.0 

ABSOLUTE ERROR FOR P-C-I, P-C-II, P-C-III, and P-C-IV 

Double Precision 



1/h 


S E 


e m 


e n 


S HE 


128 


-50.844 


-9.503- 10" 3 


-47.611 


-32.764 


64 


-107.099 


7. 761-10' 2 


-99.681 


-1.885 


32 


-214.685 


6.986 


-198.306 


18518.873 


16 


-436.640 


-12.200 


-373.796 


382.107 


8 


-898.232 


-55.310 


-660.199 


-3060.609 


4 


-1940.153 


-210.374 


-1051.118 


-2279.396 


2 


-4695.585 


-879.472 


-1475.333 


-4096.870 


CPTS 


5.65 


4.40 


5.12 


7.11 
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to h = 2 1 . The Milne method provided a valid numerical 
solution starting at h = 2 ^ and smaller, but the other 
methods still are unreliable. 

From Table 17-B, only the Hamming and Fourth Order Adams 

method showed a valid numerical solution for h 2 ^ and 

smaller, all others need much smaller h than the lowest 
. 7 

value, 2 chosen. Again on the average comparing the Milne, 
Hamming, and the Fourth Order Adams, for use in solving 
ODE VI with values of h 2 the Hamming is the first 
choice, Milne second and the Fourth Order Adams. As in 
the case of ODE II, the h values are governed by the maximum 
magnitude of C = f (x,y). If C is large then h must be 
chosen to be very small to satisfy the bounds for stability 
as had been established before. 

G. ODE VII 

True solution values are shown in Table 18 from x = 1.0 
to x = 10.0. Table 18-A showed that the Euler and Nystrom 
provide the accurate numerical solutions desired while the 
Milne and Hermite are unreliable. By analyzing the behavior 
of the roots of the Milne and Hermite methods shown in 
Table 18-B, it was observed that both exhibited unstable 
solutions as evidenced by the uniform growth of the error 
in one direction while that of the Nystrom and Euler had 
decreasing error as the range of iteration increases. Be- 
tween the Euler and Nystrom, the former produced more accu- 
rate results though it was a fraction of a second longer in 
computing time than the Nystrom. 
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TABLE 1 7 - B 



ODE VI at x = 10.0 

ABSOLUTE ERROR FOR P-C-V, P-C-VI, P-C-VII, and P-C-VIII 

Double Precision 



1/h 


1 HA 


E A2 


£ A3 


Z A4 


128 


3.848* 10' 3 


-48.786 


-16.423 


-9.371-10" 2 


64 


1.073* 10' 1 


-99.150 


-38.845 


-9.692-10" 1 


32 


4.256 


-191.895 


-48.587 


-4.337 


16 


4.556 


-347.874 


-45.839 


-34.281 


8 


-26.070 


-557.063 


188.851 


-100.571 


4 


-129.471 


-6131784 


295.289 


-255.149 


2 


-866.001 


-52.022 


3621.438 


-897.303 


OPTS 


5.0 


6.68 


7.18 


4.79 
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TABLE 18 



TRUE SOLUTION VALUES FOR ODE VII 
Double Precision 



X 


^exact 


1.0 


7 . 615941 * 10 1 


2.0 


9.640275-10" 1 


3.0 


9.950547'10‘ 1 


4.0 


9.993292-10 -1 


5.0 


9 . 999092 • 10" 1 


6.0 


9.999877'10 -1 


7.0 


9.999983-10' 1 


8.0 


9.999997-10" 1 


9.0 


9.999999-10' 1 


10.0 


9 . 999999 " 10 _1 
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TABLE 18 -A 



ODE VII at x = 10.0 

ABSOLUTE ERROR FOR P-C-I, P-C-II, P-C-III, and P-C-IV 

Double Precision 



1/h 


e e 


M 


h 


e he 


128 


5 . 91 " 10 ” 2 


-1.18551- 10 " 9 


- i ? 

-1.56-10 


-1. 750057 -10' 5 


64 


2 . 092 ' 10 11 


-6. 824836' 10" 8 


-6.78-10' 12 


-2. 691970’ lO 3 


32 


7.349-10' 11 


-5.977050-10' 6 


-3. 051- 10' 11 


-3.963316-10' 5 


16 


2.5577-10" 10 


-6 . 000472 • 10" 5 


-1.4540- 10' 10 


-1.511858-10' 5 


8 


8 . 8764- 10' 10 


-3. 322982- 10' 4 


-7.3452 • 10' 10 


-9.699544’ 10' 4 


4 


3.58311- 10" 9 


3. 314797’ 10' 3 


-3.46305-10" 9 


-6. 187207 -10' 2 


2 


-1.07753- 10' 9 


1.429417 


-5.593088-10' 7 


3. 852296- 10' 1 


CPTS 


4.42 


4.75 


4.14 


6.63 
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TABLE 18-B 



ODE at h = 0.5 

ABSOLUTE ERROR FOR P-C-I, P-C-II, P-C-III, P-C-IV 

Double Precision 



X 


S E 


Z M 


e n 


V 

^HE 


1.0 


1 . 153 ' io " 2 


- 


- 4 . 193 * 10" 5 


- 2 . 027 * 10' 4 


2.0 


8 . 578 - 10' 3 


1 . 007 - 10' 3 


- 1 . 076 * 10" 2 


- 6 . 422 - 10" 4 


3.0 


2 . 967 - 10' 3 


- 3 . 266 * 10' 3 


- 6 . 486 * 10” 3 


- 6 . 869 - 10' 4 


4.0 


7 . 349 * 10' 4 


- 1 . 193 * 10~ 2 


- 1 . 726 * 10 3 


- 2 . 874 - 10 -3 


5.0 


1 . 589 * 10" 4 


- 2 . 764 * 10' 2 


- 2 . 717 * 10' 4 


- 8 . 340 * 10' 3 


6.0 


3 . 208 - 10' 5 


- 6 . 232 - 10' 2 


- 2 . 363 * 10' 5 


- 2 . 261 * 10 ~ 2 


7.0 


6 . 221 - 10" 6 


- 1 . 250 * 10" 1 


- 1 . 228 * 10" 6 


- 5 . 534 * 10 -2 


8.0 


9 . 152 ' 10" 7 


- 1 . 906 * 10 _1 


- 1 . 132 * 10' 6 


- l . ooo - io ' 1 


9.0 


- 1 . 828 ’ 10 " 8 


- 3 . 584 * 10' 1 


- 9 . 497 * 10" 7 


- 3 . 718 - 10" 1 


10.0 


- 1 . 077* 10 " 9 


1.429 


- 5 . 593 - 10' 7 


3 . 852 * 10" 1 
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Results from Table 18-C showed that the last four methods 

produced valid numerical solutions with Hamming and Fourth 

Order Adams competing for best accuracy. As h increases 
- 7 - 5 

from h = 2 to h = 2 the Hamming method is most precise 

but from h = 2 ^ to h = 2 ^ the Fourth Order Adams method 

showed better accuracy. Computing time for the fourth Order 

Adams method is less than the Hamming method. The Euler 

method compared with these two methods as h increases, 

- 7 - 4 

started at lesser accuracy from h = 2 to h = 2 but 
maintained its accurate results up to h = 2 ^ and yielded 
a much more accurate numerical solution at h = 2 It also 

had the least computing time. Thus for classes of ODE 
belonging to ODE VII, the recommended order of choice of 
methods is as follows: Euler, Nystrom, Fourth Order Adams, 

Hamming, Third Order Adams, and lastly the Second Order 
Adams. The Milne and Hermite are not considered and should 
not be used for this particular type of ODE since they are 
both unstable. 

H. ODE VIII 

Table 19 presents the true solution values for ODE VIII 
from x = 1.0 to x = 10.0. Data analysis from Tables 19-A 
and 19-B indicates that any of the eight P-C sets can be 
used to solve ODE VIII if no specific criterion for accuracy 
is needed, since each method yields- good and valid numerical 
solutions. However, for purposes of comparison, the Hamming 
method clearly stands out to be the best in terms of accuracy 
with Milne coming in second. The Milne method provides the 
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TABLE 18- C 



ODE at x = 10.0 

ABSOLUTE ERROR FOR P-C-V, P-C-VI, P-C-VII, and P-C-VIII 

Double Precision 



1/h 


Z HA 


E A2 


M 

> 


^A4 


128 


0 


- 2 . 398169 * 10" 7 


-5.968010’ 10" 10 


0 


64 


0 


-1.929706-10" 6 


-1.07296"10" 9 


0 


32 


0 


- 1 . 56 1650 * 10 " 8 


-1.87211* 10 " 9 


-5.0- 10' 14 


16 


5 .665862- 10" 8 


-1.277970- 10' 4 


-3.04922-10' 9 


- 1 2 

-1.36-10 


8 


5 . 374510- 10' 8 


-1.068471-10" 3 


-1.119715-10' 8 


1 . 715985 • 10" 7 


4 


2.632371-10' 7 


-9. 306166’ 10' 3 


-3.729009’10' 7 


-1.727044-10' 7 


2 


1. 049942- 10' 2 


-8. 771719* 10 -2 


-2. 033756- 10' 3 


-9.570674-10' 3 


OPTS 


5.15 


7.05 


6.20 


4.82 
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TABLE 19 



TRUE SOLUTION VALUES FOR ODE VIII 
Double Precision 



X 


^exact 


1.0 


2.319776 


2.0 


2.482577 


3.0 


1.151562 


4.0 


4 . 691641 - 10' 1 


5.0 


3 . 833049 * 10' 1 


6.0 


7 . 562256 * 10 _1 


7.0 


1.928970 


O 

CO 


2.689507 


9.0 


1.510013 


10.0 


5 . 804096 * 10 _1 
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TABLE 19 -A 



ODE at x = 10.0 

ABSOLUTE ERROR FOR P-C-I, P-C-II, P-C-III, and P-C-VI 

Double Precision 



1/h 


Z E 


e m 


e n 


Z HE 


128 


1. 391257* 10' 6 


-9 . 332 * 10 " 11 


7.611170* 10" 7 


1. 595845* 10' 7 


64 


6. 532096* 10' 6 


-2.17509-10" 9 


5.204737-10" 6 


5.639547*10" 7 


32 


1 . 850381 * 10” 5 


-5. 498490* 10" 8 


1.836626* 10‘ 5 


4.499001* 10" 6 


16 


7. 305628* 10" 5 


-6 . 440709 * 10~ 7 


7.254094* 10' 5 


3. 580536 * 10” 5 


8 


2.827282-10' 4 


-4.145628-10" 6 


2. 733740* 10" 4 


2 . 841152 * 10" 4 


4 


9. 532272* 10" 4 


-9. 793708* 10' 5 


7 . 115245 * lO 4 


2. 270258-10" 3 


2 


1 .051152* 10" 3 


-2 . 707671* 10' 3 


-1.087262*10" 2 


1. 866206* 10" 2 


OPTS 


6.14 


4.63 


5.11 


6.48 
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TABLE 19 -B 



ODE at x = 10.0 

ABSOLUTE ERROR FOR P-C-V, P-C-VI, P-C-VII, and P-C-VIII 

Double Precision 



1/h 


Z HA 


Z A2 


Z A3 


S A4 


128 


- 1 3 

-9.7*10 ^ 


9.902758-10" 7 


2 . 367453 * 10" 7 


-2.8983-10" 10 


64 


-4.057-10" 11 


3.303524-10" 6 


1. 886624* 10' 6 


-5.53331*10" 9 


32 


-1. 838* 10" 9 


7.843814-10' 6 


1. 498027- 10' 5 


-1.155997-10' 7 


16 


-1.006735* 10" 7 


-1.359301-10' 5 


1.181971* lO -4 


-1.563035-10' 6 


8 


-1.414250-10" 6 


-4.482144* 10' 4 


9.248289* 10" 4 


-2. 197783* 10" 5 


4 


-7.603359* 10" 5 


-6.125271-10' 2 


5. 783513* 10" 2 


5. 113170- 10' 3 


2 


-9.956171’10 -4 


-6.125271*10' 2 


5 . 783513* 10" 2 


5.113170’10" 3 


OPTS 


5.31 


7.61 


8 . 11 


5.09 
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least computing time with Hamming second. Thus for problems 
of the type ODE VIII, the Hamming method is the most highly 
recommended, followed by the Milne, Fourth Order Adams, 

Euler, Hermite, Nystrom, Third Order Adams, and lastly the 
Second Order Adams. 

I. ODE IX 

True solution values are listed in Table 20 for ODE IX 
for the range of integration x = 1.0 to x = 10.0. From 
Tables 20-A and 20-B it was observed by analyzing the 
results that all eight methods provided valid numerical solu- 
tions. In terms of better accuracy on the average and the 
least computing time the Milne method is the best candidate 
with Hamming coming in next, followed by the Fourth Order 
Adams, Hermite, Third Order Adams, Euler, Nystrom and 
lastly the Second Order Adams. Thus if choice is to be 
made for solving classes of ODE belonging to ODE IX, the 
Milne method is highly recommended with Hamming as the next 
alternate . 

J. ODE X 

This differential equation is an example of a controlled 
variable problem, where one variable is expressed as a func- 
tion of x and substituted into the differential equation of 



the other variable to 
original two variable 


form a single first order ODE. 
equations are: 


The 


y{ = i.38 yi - 


0 . 8ly 2 


(2-74) 


7 2 = 2 *i6y 1 - 


1.92y 2 


(2-75) 



142 



TABLE 20 



TRUE SOLUTION VALUES FOR ODE IX 
Double Precision 



X 


^exact 


1.0 


2.069535 


2.0 


2.249705 


3.0 


4.179309 


4.0 


9.462527 


5.0 


10.633343 


6.0 


17.564095 


7.0 


42.421352 


8.0 


50.806493 


9.0 


74.608406 


10.0 


186.463649 
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TABLE 2 0 -A 



ODE IX at x = 10.0 

ABSOLUTE ERROR FOR P-C-I, P-C-II, P-C-III, and P-C-IV 

Double Precision 



1/h 


Z E 


M 


Z N 


e he 


128 


-1.591-10" 3 


1. 791* 10" 8 


-1 .658-10" 3 


1.403-10" 5 


64 


-6 . 435 ‘ 10" 3 


1.666-10" 7 


-6.522* 10" 3 


1.173-10" 4 


32 


-2. 578-10' 2 


2 . 341 * 10" 6 


-2.636-10' 2 


9.226-10" 4 


16 


-1.025"10 _1 


3.806 - 10" 5 


-1.079- 10 _1 


7.170-10' 3 


8 


-4.022'1- -1 


1.016-10' 3 


-4.569-10' 1 


5 . 390 * 10" 2 


4 


-1.450 


4 . 234 ' 10" 2 


-2.090 


3 . 761 ’ 10 " 1 


2 


-3. 737 


9. 791* 10" 1 


-10.452 


2 .559 


OPTS 


6.99 


5.40 


7.09 


6.42 
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TABLE 20 -B 



ODE IX at x = 10.0 

ABSOLUTE ERROR FOR P-C-V, P-C-VI, P-C-VII, and P-C-VIII 

Double Precision 



1/h 


S HA 


£ A2 


V* 

^A3 


E A4 


128 


-2. 795*10" 9 


-1.648* 10" 3 


6. 166* 10" 5 


6 . 922 * 10" 8 


64 


-6 . 428 * 10" 8 


-6.720-10' 3 


4 . 827 * 10 4 


9. 598* 10" 7 


32 


-2.136-10" 6 


-2. 788* 10~ 2 


3. 693 * 10 ” 3 


1.610* 10' 5 


16 


-8.656-10' 5 


-1.190-10 -1 


2.690-10' 2 


2.948-10' 4 


8 


-9 . 897 * 10~ 4 


-5. 288* 10" 1 


1. 755* 10' 1 


7 . 000 * 10" 3 


4 


-5.696-10" 2 


-2 . 433 


8.959*10 _1 


2 . 261 " 10 1 


2 


2.108 


-10.866 


3. 597 


4.380 


CPTS 


5. 73 


7.93 


8.33 


5.44 
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with initial conditions 



y 1 (0) = -2.9995 



(2-76) 



y 2 (0) = 4.0010 



(2-77) 



and the analytic solutions given by 



y^ = 0.0005 e 



- 3x 



3e 



- 0 . 3x 



(2-78) 



y 2 = 0.001 e 



- 3x 



+ 4e 



-0 . 3x 



(2-79) 



Since the analytical solutions of both differential equations 
are known, it is possible to treat one variable as a function 
of x and substitute it to the differential equation of the 
other variable to form a single first order ODE which now can 
be solved by the P-C sets. Choosing y 2 as the known function 
given by (2-79) and substituting it in (2-74), the resulting 
equation is ODE X with initial condition given by (2-76) and 
exact solution by (2-78). Note that y 2 remains a function 
of the single independent variable x and as such changes as 
x changes. The other variable y^ forms a first order dif- 
ferential equation in the standard form y| = f(x,y). The 
point that is being illustrated here is that this concept can 
be applied to problems like rate of chemical reaction, falling 
bodies, aircraft or ballistic missiles flight wherein time 
variable is the most important consideration. All other 
variables can be expressed as functions of the single inde- 
pendent variable time and given constants or initial condi- 
tions. In such cases the problem can be reduced to a single 
differential equation and the methods considered herein can 
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be applied (time (t) variable as x) . This concept does not 
discount the fact that these P-C methods can be extended to 
solve simultaneous differential equations with some modifi- 
cations in the algorithms. However, since this extension 
is not relevant to the purpose of this particular paper in 
analyzing and comparing the predictor-corrector methods as 
applied to a variety of ODEs , it is not considered here, 
but will be mentioned as a further field of study. 

Returning now to the solution of ODE X, Table 21 shows 
the true solution values from x = 1.0 to x = 10.0. From 
Table 21-A the Milne seemed to produce the best accuracy 
but by analyzing further the behavior of these errors in 
step-by-step integration, which is shown in Table 21-B, it 

i 

is noted that the errors of the Milne method are increasing 
in a uniform fashion such that as the range of integration 
is increased the error grows while the true solution values 
as shown in Table 21 decrease. Thus eventually the error 
will overcome the solution. The same is true for the 
Hermite method, while for the Euler and the Nystrom the 
error decreases as the range of integration increases thus 
providing a valid numerical solution. The behavior of the 
Milne and Hermite methods is expected since for ODE X it is 
easy to see that C < 0, and as such the Milne and Hermite 
methods are unstable. Between the Euler and the Nystrom, 
the Euler method is an easy choice both in accuracy and 
computing time. 
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TABLE 21 



TRUE SOLUTION VALUES FOR ODE X 
Double Precision 



X 


^exact 


1.0 


-2.222429 


2.0 


-1.646433 


3.0 


-1.219708 


4.0 


-9.035826-10" 1 


5.0 


-6.693904-10' 1 


O 

VO 


-4.958966-10' 1 


7.0 


-3.673692 * 10 _1 


8.0 


-2. 721538’ 10' 1 


9.0 


-2 . 016165 * 10* 1 


10.0 


-1.493612* 10' 1 
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TABLE 21 -A 



ODE at x = 10.0 

ABSOLUTE ERROR FOR P-C-I, P-C-II, P-C-III, and P-C-IV 

Double Precision 



1/h 


Z E 


E M 


e n 


e he 


128 


4.839775-10' 7 


-1.20S-10' 11 


-3.684268-10" 8 

\ 


-9.818242-10' 4 


64 


-9. 887212 * 10" 8 


-5. 3459*10" 10 


-8. 584152* 10' 8 


-3.913399-10' 3 


32 


-4.893782-10' 7 


-3.606138-10" 8 


-3.832304*10" 7 


-1.558160-10" 2 


16 


-1.149738-10' 6 


-4. 286345* 10~ 7 


-1. 793736* 10' 6 


-6. 236379* 10" 2 


8 


-3. 757046*10 _6 


-1. 5761S3* 10' 5 


-4.903244-10' 6 


-2.594684-10' 1 


4 


-9.800289* 10~ 8 


-2. 583136- 10" 4 


-2.093989-10" 5 


-1.288607 


2 


2 . 784363 * 10" 4 


- 8 . 4 36503 * 10 ~ 3 


-1. 350331 * 10" 4 


-11.095378 


CPI'S 


7.03 


5.33 


10.0 


9.83 
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TABLE 21- B 



ODE X with h = 0.5 

ABSOLUTE ERROR FOR P-C-I, P-C-II, P-C-III, and P-C-IV 

Double Precision 



X 


Z E 






2 he 


1.0 


2.187-10' 3 




-5.056-10" 4 


-5.661-10" 2 


2.0 


2.565-10" 3 


-3.821-10' 5 


-1. ioi-io' 3 


-8.576-10' 2 


3.0 


2.144-10' 3 


- 1 . 905 * lO 4 


-1. 016 " 10" 3 


-1.542-10' 1 


4.0 


1 . 651 ' 10" 3 


-3. 765- 10" 4 


-7.984-10' 4 


-2.833-10" 1 


5.0 


1. 239-10" 3 


-6.420- 10" 4 


-6.013-10" 4 


-5. 219- 10' 1 


6.0 


9. 223* 10" 4 


-1.076* 10 " 3 


-4.475*10 _4 


-g.oio'io' 1 


7.0 


6 . 843" 10" 4 


-1. 800* 10' 3 


-3.319-10' 3 


-1.772 


8.0 


5.072- 10" 4 


-3.012- 10' 3 


-2.460-10' 3 


-3.266 


9.0 


3. 758* 10' 4 


-5.041*10 3 


-1.822-10' 3 


-6.020 


10.0 


2 . 784 ’ 10 ~ 4 


- 8 . 436 ' 10 ” 3 


-1. 350- 10' 3 


-11.095 
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- 




From Table 21-C, all the last four methods exhibited 
stable numerical solutions. For accuracy the Hamming method 
offers the best choice, with Fourth Order Adams coming 
next, though the computing time for Hamming is a little 
longer than the Fourth Order Adams. Therefore for problems 
in the class of ODE X, the order of recommended preference 
is: Hamming, Fourth Order Adams, Third Order Adams, Euler, 

Nystrom, and the Second Order Adams as the last choice. The 
Milne and Hermite methods should not be used for this par- 
ticular type of problem. 

K . SUMMARY 

- 7 - 1 

For series of h values from h = 2 to h = 2 and range 
of integration up to x = 10.0, Tables 22, 22-A, and 22-B 
list the summary of results for the eight predictor-corrector 
sets considered, each run on the test ODEs I to X. From 
these comparative results, it is clearly established that 
the best numerical methods in order of decreasing efficiency 
in general are: 

1. Hamming Modified P-C Set 

2. Fourth Order Adams -Moulton P-C Set 

3. Milne P-C Set 

4. Third Order Adams -Moulton P-C Set 

5. Euler P-C Set 

6. Nystrom P-C Set 

7. Second Order Adams -Moulton P-C Set 

8. Hermite P-C Set. 
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TABLE 21- C 



ODE X at x = 10.0 

ABSOLUTE ERROR FOR P-C-V, P-C-VI, P-C-VII, and P-C-VIII 

Double Precision 



1/h 


y 

l ha 


E A2 


Z A3 


Z A4 


128 


0 


-1 . 149665 * 10" 7 


-9.07653* 10" 9 


-4.0-10' 14 


64 


0 


-8. 476336- 10' 7 


-7. 319508-10' 8 


- 13 

-6.9-10 10 


32 


9 . 0 ' 1 0 ~ 14 


-6. 538981- 10' 6 


-5 . 949594 * 10" 7 


-1 . 380 * 10" 11 


16 


8. 773099 -10' 8 


-5. 207967-10' 5 


-4.912304-10" 6 


-3.0945-10" 10 


8 


-7.06S21S-10' 8 


-4. 272440* 10' 4 


-4. 180375-10' 5 


-2.697520-10* 8 


4 


4.17081-10' 9 


-3.630456* 10" 3 


-3. 757718’10' 4 


1 . 5 5 361 4 * 1 0 ~ 7 


2 


-4. 137219- 10" 5 


-3. 165651- 10' 2 


- 3. 628554 ’ 10" 3 


-1 . 293084 * 10 " 6 


CPTS 


6.02 


9.45 


9.98 


5.45 
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TABLE 22 



SUMMARY OF RESULTS FOR ODE I TO ODE IV 
FOR h= 2 ~ 7 TO h=2 ~ 1 AT x=10.0 



ODE 

Number 


Stable P-C Sets in 
Order of Preference 


Average 

Accuracy 


Computing 
Time 
in sec. 


Unstable 
P-C Sets 


I 


1. 


Hamming 


IO' 4 


5.76 


1. Milne 




2. 


Fourth Order Adams 


10~ 3 


5.99 


2. Hermite 




3. 


Third Order Adams 


10" 3 


8.55 






4. 


Euler 


10- 2 


7.65 






5. 


Nystrom 


10- 2 


7.46 






6. 


Second Order Adams 


10' 2 


7.98 




II 


1. 


Hamming 


10- 4 


5.26 


1. Hermite 




2. 


Fourth Order Adams 


10' 3 


4.96 






3. 


Milne 


10' 3 


4.69 






4. 


Euler 


10 _1 


6.38 






5. 


Second Order Adams 


10' 1 


7.79 






6. 


Third Order Adams 


10' 1 


8.36 






7. 


Nystrom 


IO -1 


5.34 




III 


1. 


Hamming 


10~ 8 


6.23 


None 




2. 


Milne 


10' 7 


5.59 






3. 


Fourth Order Adams 


10' 7 


8.23 






4. 


Third Order Adams 


10 _5 


10.0 






5. 


Hermite 


10' 4 


7.92 






6. 


Euler 


IO -4 


7.81 






7. 


Second Order Adams 


10~ 4 


10.0 






8. 


Nystrom 


10" 4 


7.12 




IV 


1. 


Hamming 


IO' 10 


4.45 


None 




2. 


Mi lne 


10' 9 


3.83 






3. 


Fourth Order Adams 


10' 8 


4.02 






4. 


Third Order Adams 


10' 6 


7.02 






5. 


Nystrom 


io" 5 


3.67 






6. 


Second Order Adams 


10- 5 


6.46 






7. 


Euler 


io' 5 


4.39 






8. 


Hermite 


10' 4 


6.74 
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TABLE 22 -A 



SUMMARY OF RESULTS FOR ODE V TO ODE VIII 
FOR h= 2 " 7 TO h=2 _1 AT x=10.0 



ODE 


Stable P-C Sets in 




Computing 


1 

Unstable 


Number 


Order of Preference 


Accuracy 


i line 
in sec. 


P-C Sets 


V 


1 . Hamming 


lO' 13 


6.31 


None 




2. Milne 


10- 12 


5.58 






3. Fourth Order Adams 


10- 12 


7 .78 






4. Third Order Adams 


10' 11 


10.0 






5. Hermite 


icf 10 


5.35 






6. Ny strom 


10' 7 


5.16 






7. Euler 


lO' 7 


5.63 






8. Second Order Adams 


10' 7 


10.0 




VI 


1. Hamming 


10~ 2 


5.0 


1. Euler 




2. Fourth Order Adams 


10' 1 


4.79 


2. Nystrom 




3. Milne 


10' 1 


4.40 


3. Hermite 










4. Third Order 










Adams 










5. Second Order 










Adams 


VII 


1. Euler 


l(f 10 


4.42 


1. Milne 




2. Ny strom 


10~ 9 


4.14 


2. Hermite 




3. Fourth Order Adams 


10" 8 


4.82 






4. Hamming 


10" 8 


5.15 






5. Third Order Adams 


io" 7 


6.20 






6. Second Order Adams 


10' 4 


7.05 




VIII 


1. Hamming 


io' 8 


5.31 


None 




2. Milne 


10' 7 


4.63 






3. Fourth Order Adams 


io' 6 


5.09 






4. Euler 


io” 5 


6.14 






5. Hermite 


io' 5 


6.48 






6. Nystrom 


io' 4 


5.11 






7. Third Order Adams 


io" 4 


8.11 






8. Second Order Adams 


10‘ 4 


7.61 
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TABLE 22-B 



SUMMARY OF RESULTS FOR ODE IX AND ODE X 
FOR h=2 " 7 TO h=2 _1 AT x=10.0 



ODE 

Number 


Stable P-C Sets in 
Order of Preference 


Average 

Accuracy 


Computing 

Time 

in secs. 


Unstable 
P-C Sets 


IX 


1. 


Milne 


IQ' 4 


5.40 






2. 


Hamming 


10* 4 


5.73 






3. 


Fourth Order Adams 


io ' 4 


5.44 






4. 


Hermite 


10" 3 


6.42 






5. 


Third Order Adams 


10- 3 


8.33 






6. 


Euler 


io " 2 


6.99 






7. 


Nys trom 


10- 2 


7.09 






8. 


Second Order Adams 


io ' 2 


7.93 




X 


1. 


Hamming 


icf 10 


6.02 


1. Milne 




2. 


Fourth Order Adams 


10' 9 


5.45 


2. Hermite 




3. 


Third Order Adams 


10- 6 


9.98 






4. 


Euler 


10" 6 


7.03 






5. 


Nys trom 


IO' 6 


10.0 






6. 


Second Order Adams 


io - 5 


9.45 
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The most important criteria applied is accuracy followed by 
computing time. It is obvious that the Hamming method did 
not produce the least computing time as compared to the 
other methods especially with that of Milne and the Fourth 
Order Adams methods. This difference of less than a second 
in computing time is attributed to the fact that the Hamming 
method used additional computation for the truncation error 
to modify the predicted and corrected values for every 
iteration. Analyzing the computing time for each method 
within the range of their stability, the list below shows 
the order of increasing computing time: 





P-C Set 


Average 

Computing Time 


1 . 


Mi lne 


4.82 


secs . 


2. 


Hamming 


5.60 


secs • 


3. 


Nys trom 


5.65 


secs . 


4. 


Fourth Order Adams 


6.11 


secs . 


5. 


Euler 


6.19 


secs • 


6. 


Hermite 


6.58 


secs • 


7. 


Second Order Adams 


8.40 


secs . 


8. 


Third Order Adams 


8.69 


secs . 



This confirmed the idea that convergence of a corrector does 

not necessarily assure the convergence towards the true 

solution values, y(x ), but only to some definite value, 

y . , . This is best illustrated in the case of Hermite and 
' n+1 

Third Order Adams methods with 6.58 secs and 8.69 secs average 
computing time, respectively. By recalling from the order 
of efficiency of the P-C sets, the Third Order Adams ranks 
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fourth in accuracy while Hermite ranks eighth, which points 

to the fact that though the Hermite method converges much 

faster than the third order Adams, its accuracy is much 

poorer than the Third Order Adams. Clearly then it can be 

seen that the Hermite method is converging only to some 

definite value, y + ^, but not to the true solution, y(x n ), 

otherwise it should have much, better accuracy than the Third 

Order Adams method, which took, much longer to converge to 

a definite value y in each iteration. This analysis is 

meaningful due to the fact that a standard mode of applica- 

2 

tion is used, that of P(EC) . It should be noted however 
that rapid convergence is essential to accuracy as shown 
by the Hamming, Milne, and Fourth Order Adams methods which 
all rank well above the other methods in both accuracy and 
computing time in general. 
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X. CONCLUSIONS 



The following conclusions can be drawn and recommenda- 
tions made from the analysis and numerical results obtained 
in this paper. 

1. Numerical instability results from extraneous solutions 

of the difference equations which bear no connection to 

the exact solution. The conditions of asymptotic (strong) 

and absolute stability can be seen clearly by reference 

to the extraneous solutions. If in the limit h -> 0 , the 

extraneous solutions vanish as n 00 , then the method is 

strongly stable and convergent. If, for values of h less 

than some h , the extraneous solutions vanish as n -> <», the 
o 

method is absolutely stable. Relative stability is a sig- 
nificant concept for ODE where f >0. The condition pro- 
vides that extraneous solutions will not grow more rapidly 
or decay more slowly than the true solution. Thus, before 
a method is used, the characteristic roots of its related 
difference equation must be analyzed. 

2. The finite real negative stability bounds for the P-C 
methods considered, as determined by numerical experiments, 
are in reasonable agreement with theoretical predictions. 

As such experimental bounds could be used as a guide to the 
proper selection of a method to solve a particular problem. 

3. The stability of a P-C set depends on both the predictor 
and corrector equations. Thus, when two P-C sets have the 



same corrector but different predictors, choose the one with 
a better predictor. Likewise, when two P-C sets use the 
same predictor with different correctors, choose the one 
with a better corrector. 

4. P-C sets with higher order truncation errors are not 
necessarily better. In choosing a P-C set, the order of 
truncation error should not be the sole criterion. 

5. The convergence of the corrector formula will ensure 
that the sequence of approximations will converge to some 
definite value but not necessarily to the true solution 
values. As such, the fast computation time of a method 
does not necessarily imply greater accuracy for the result- 
ing numerical solution. 

6. If the integration will involve a large number of steps, 
a stable method should be used. 

7. If the function evaluation is lengthy, the range of inte- 
gration large, and better accuracy and lesser cost of com- 
puting time are prime considerations, then the following P-C 
methods are recommended, based on the overall efficiency 
they have demonstrated in numerical experiments on a wide 
variety of different test ODEs : 

a) Hamming P-C Sets 

b) Fourth Order Adams P-C Set 

c) Milne P-C Set 

d) Third Order Adams P-C Set 

e) Euler P-C Set 

f) Nystrom P-C Set 
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g) Second Order Adams P-C Set 

h) Hermite P-C Set 

The listing indicates order of preference. However caution 
must be exercised not to use Milne and Hermite methods if 
f <0. 

y 

8. For maximum accuracy to be achieved, select a step size 
h, for which the truncation error and the roundoff error have 
equal orders of magnitude. 

9. If high accuracy results with less roundoff error are 
desired, the procedure should use double precision, which 
involves only slightly more computing time than single pre- 
cision in an IBM 360/67 computer. 
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XI. EXTENSIONS 



The following topics, which are direct extensions of 
this paper, are interesting areas of further research and 
study in the wide field of numerical analysis: 

1. Stability versus Accuracy 

A necessary condition for numerical integration is 
stability' that is, a stable value of h must be employed 
when a fnite stability boundary exists. Nevertheless, 
there is no guarantee that all values of h from zero to the 
limiting \aLue will yield accurate results. Thus the 

I 

question that must be faced is: Of what value is a method 

that is stable in a region where it is inaccurate? The 
answer to this question will show that the relationship 
between stability and accuracy is fundamental to the choice 
of h for a particular method. 

2. Increasing the Stability Bounds 

It has been shown that in some methods the real stability 
bound is somewhat constrained. Two possible ways where the 
stability bounds can be increased are: 

a) By the use of aweighting factor which involves experi- 
menting with different values of the free parameters that 
are used in the method by undetermined coefficients, similar 
to Hamming's corrector formula derivation. 

b) By averaging, which means computing the new value of 
y n+1 after a finite number of steps (say fifty) of a particular 
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method, as the average of the old y n+ ^ and a value called 
y n+ i generated by another method. This should be done 
periodically . 

3. Predictor-Corrector Methods Interaction 

The interaction of the different methods in solving a 
particular problem involves use of the P-C sets as sub- 
routines, wherein one method will be used to provide the 
solution for smaller values of h, then another method will 
be used to solve for larger values of h. This should be 
quite interesting since it has been shown that some methods 
exhibit better accuracy for small values of h, then become 
inaccurate for large values of h, while other methods main- 
tain their accuracy up to large values of h but might be 
prohibitive to use for small values of h due to their longer 
computing time. 

4. Solution of Systems of Differential Equations 
Revision of the algorithms to enable the solution of 

systems of differential equations and higher order differ- 
ential equations. This extension is straightforward. It 
involves only extra computational steps using the same 
formulas. Complexity arises from the need to use vectors 
and arrays for temporary storage. 

5. Practical Applications 

Applications to real-life problems such as heat flow 
problems, simple electrical circuits, force problems, rate 
of bacterial growth, rate of decomposition of radioactive 
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material, crystallization rate of a chemical compound, rate 
of population growth, and so forth. 

6. Adjustment of the Step Size During the P-C Solution 

This involves monitoring the modified truncation error 
developed by Hamming. If the value is too large, then 
the step size is too large, and the calculation should be 
repeated with a smaller value of h, say h/2. Note that it 
is not necessary to go back to the beginning of the calcula- 
tion but only to the point at which the truncation error 
becomes too large. If the truncation error is too small, 
computation time is being wasted and h should be increased, 
say to 2h. 
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FLOW CHART FOR ALGORITHM 1 
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FLOW CHART FOR ALGORITHM 2 
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PROGRAM ONE 



EULER PREDICTOR-CORRECTOR METHOD. ODE: Y ' =-Y 
WITH Y ( 0 ) =1 . TO SOLVE ANOTHER ODE SIMPLY CHANGE 
INITIAL CONDITIONS AND FUNCTION SUBROUTINES. 



C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



c 

c 

c 

5 

c 

c 

c 



c 

c 

c 

c 

c 

c 

7 

c 

c 

c 






* FORTRAN PROGRAM FOR EULER PREDICTOR-CORRECTOR 

* METHOD IN THE NUMERICAL SOLUTION OF THE ORDINARY 

* DIFFERENTIAL EQUATION IN THE FORM OF DY/DX=F(X,Y) 

* WITH THE INITIAL CONDITION Y(XO)=YO 

* THE PARAMETERS TO THE PROGRAM HAVE THE FOLLOWING 



* MEANINGS... 

* XMAX THE TERMINAL BOUNDARY CONDITION 

* EPS THE CONVERGENCE TEST CONSTANT 

* MAX ITHE MAXIMUM NUMBER OF ITERATIONS 

* YP_I THE PREDICTED VALUE OF THE NEXT 

* SOLUTION POINT 

* YC THE CORRECTED VALUE OF THE NEXT 

* SOLUTION POINT 

* FP THE PREDICTED FUNCTION EVALUATION OF YP 

* FC THE CORRECTED FUNCTION EVALUATION OF YC 

* H ITHE STEP SIZE TO BE USED TO ADVANCE 

* THE POINT OF INTEGRATION 

* ICON THE DESIRED POINT OF PRINTING THE 

* COMPUTED SOLUTION FOR EVERY FIXED 

* NUMBER OF INTEGRATION POINTS 

* YEXACT THE TRUE SOLUTION VALUES 

* ERROR THE ABSOLUTE DEVIATION OF THE NUMERICAL 

* SOLUTION FROM THE TRUE SOLUTION 

* FCT FUNCTION SUBROUTINE TO EVALUATE ODE 

* FUNCTION 

* EXACT FUNCTION SUBROUTINE TO COMPUTE TRUE 

* SOLUTION VALUES 



sV 4. V. 'V .V O'- J/ O. U. 4/ -JU 4. 4. - 

^ ^ ^ ^ •y*' -y y .r. * **| ' y c 



* *r* 'r v 'r 



4. 4> 4. 4^ 4. 4> 4. .V 4. 4. X X 4. .V 4U 4> 4, 4. .V 4r 4U 4# 4. 



INITIALIZE INPUT CONSTANTS 

MAX= 2 
XMAX= 1 0.0 
EPS=0. 000001 

READ STEP SIZE AND CONTROL PRINTING 

READ ( 5 » 100 ) H t I CON 

SPECIFY INITIAL CONDITION 

X0=0.0 
Y0=1 .0 

TEST FOR END OF INPUT DATA 

I F { I CON— 0 ) 45,45,7 

PRINT INPUT DATA SET 

WRITEI6, 1000) H , I CON 
WR I T E ( 6 , 1 001 ) XO,YO 

COMPUTE NEXT SOLUTION POINT 

NC=( XMAX+H/2. I /H 
FP=FCT (XO, YOI 
DO 40 N=1 , NC 
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-H- -5M(- tt-M- # -K- ■* -if- -M- -if- 



non non ooo non ooo ooo non 



11 



M=1 

XO=XO+H 
YP=Y0+H*FP 
FC = FCT ( XO » YP ) 

YC=Y0+H/2.*(FP+FC) 

DELY=ABS ( YC-YP ) 

TEST FOR CONVERGENCE 

IF(DELY-EPS) 30,30,15 

TEST IF MAXIMUM NUMBER OF ITERATIONS IS SATISFIED 

15 IF(M^MAX) 20,20,30 
20 YP=Y C 

M=M+ 1 
GO TO 11 

COMPUTE TRUE SOLUTION VALUE 
30 YEXACT=EXACT (X0) 

COMPUTE ERROR IN CALCULATION 
ERR0R= YEXACT— Y C 

TEST IF DESIRED POINT OF PRINTING IS REACHED 

I F ( { N/ 1 CON* ICON) .EQ.N) WRI T E {6 , 1 002 ) XO ,YP , YC , YEXACT , 

1 ERROR 

UPDATE REQUIRED PARAMETER VALUES 

Y0=Y C 
FP=FC 

40 CONTINUE 

LOOP BACK TO READ NEXT. SET OF INPUT DATA 

GO TO 5 
45 STOP 

100 FORMAT (F10.7,I10) 

1000 FORMAT (/// /41X » • *1 MPUT DATA SET USED** //43X, 'H =', 
1F10. 7/43X , • ICON =',I10) 

1001 FORMAT I//29X, '**RESULTS OF EULER PREDICTOR-CORRECTOR 

1METHOD** ' // 10X, ' X ',5X,* YP ’,5X,' YC ', 

15X , ' YEXACT • , 5X , ' ERROR • // 12X , F10 . 7 ,44X , 
1F10.7) 

1002 FORMAT ( • 0 ' , 9X , F15 .7 , 2X , F 15 . 7 , 2X , F 1 5 . 7, 2X , F 1 5 . 7 , 2X , 

IF 1 5. 7 ) 

END 



FUNCTIGN FCT (XN,YN) 

FCT=-YN 

RETURN 

END 



FUNCTION EXACT(XX) 
EXACT=EXP( -XX ) 
RETURN 
END 
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❖INPUT DATA SET USED* 

H = 0.5000000 
ICON = 2 

❖❖RESULTS OF EULER PREDICTOR-CORRECTOR METHOD** 



X 


YP 


YC 


YEXACT 


ERROR 


0. 0000000 






1 .0000000 




1.0000000 


0.3588867 


0.3634033 


0.3678794 


0. 0044761 


2. 0000000 


C. 1309212 


0.1325720 


0.1353353 


0.0027633 


3.0000000 


0.0477613 


0.0483635 


0.0497871 


0.0014236 


4. 0000000 


0.0174238 


0. 0176435 


0.0183156 


0.0006722 


5.0000000 


0.0063564 


0.0064365 


0.0067379 


0.0003014 


6.0000000 


0.0023189 


0.0023481 


0.0024788 


0. OdOl 307 


7.0000000 


0. 0008459 


0.0008566 


0.0009119 


0.0000553 


8.0000000 


0.0003086 


0.0003125 


0.0003355 


0.0000230 


9.0000000 


0.0001126 


0.0001140 


0.0001234 


0.0000094 


10.0000000 


0.0000411 


0.0000416 


0.0000454 


0.0000038 



* I NPUT DATA SET USED* 

H = 1.0000000 
ICON = 1 

**RE SULT S GF EULER PREDICTOR-CORRECTOR METHOD** 



X 


YP 


YC 


YEXACT 


ERROR 


0.0000000 






1 .0000000 




1.0000000 


0.2500000 


0.3750000 


0.3678794 


-0.0071206 


2.0000000 


0.1562500 


0.1718750 


0.1353353 


-0.0365397 


3.0000000 


0.0507813 


0.0683594 


0.0497871 


-0.0185723 


4.0000000 


0.0253789 


0. 0300293 


0.0183156 


-0.0117137 


5.0000000 


0.0095825 


0.0122986 


0.0067379 


-0.0055606 


6.0000000 


0.0044327 


0.0052910 


0 .0024788 


-0.0028122 


7.0000000 


0.0017519 


0.0021987 


0.0009119 


-0.0012868 


8.0000000 


0.0007731 


0.0009362 


0.0003355 


-0.0006007 


9.0000000 


0.0003156 


0.0003919 • 


0.0001234 


-0.0002685 


10.0000000 


0. 0001361 


0.0001660 


0.0000454 


-0.0001206 
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PROGRAM TWO 



MILNE PREDICTOR-CORRECTOR METHOD. ODE: Y*=-Y 
WITH Y( 0 ) =1 . TO SOLVE ANOTHER ODE SIMPLY CHANGE 
INITIAL CONDITIONS AND FUNCTION SUBROUTINES. 



C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



c 

c 

c 

5 

c 

c 

c 



c 

c 

c 

c 

c 

c 

7 

c 

c 



* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

if 

* 

* 

* 

* 

if 

* 

if 

4U 

nr 

* 

* 

* 

* 

* 



4 J/ *V >*' O. <JL> 4. 4« *V 4* '** 4/ V; »<* 4* -.V «JU Jf >V J* 4. J# 4. 4. J» 4> 4f 4> 4* 4* 4« 4» * 4. 4r' 4. 4^ 

- 'V n* 'r -r ^ ^ -i- r* *r v *r *r v i' -r *>' »r J r «r» 'p ’t‘ 'r v n- -r* rr nr nr n'* nT n-- -r -V* n~ -r n' n* n' v * r *v nr nr nr *v 

FORTRAN PROGRAM FOR MILNE PREDICTOR-CORRECTOR 
METHOD IN THE NUMERICAL SOLUTION OF THE ORDINARY 
DIFFERENTIAL EQUATION IN THE FORM OF DY/DX=F(X,Y) 
WITH THE INITIAL CONDITION Y(XO)=YO 
THE PARAMETERS TG THE PROGRAM HAVE THE FOLLOWING 
MEANINGS.. . 

_THE TERMINAL BOUNDARY CONDITION 
_ THE CONVERGENCE TEST CONSTANT 
_ITHE MAXIMUM NUMBER OF ITERATIONS 

THE PREDICTED VALUE OF THE NEXT 

SOLUTION POINT 

THE CORRECTED VALUE OF THE NEXT 

50LUT I ON POINT 
THE FUNCTION EVALUATION 
“THE STEP SIZE TO BE USED TO ADVANCE 
THE POINT OF INTEGRATION 

THE DESIRED POINT OF PRINTING THE 

COMPUTED SOLUTION FOR EVERY FIXED 
NUMBER OF INTEGRATION POINTS 
YEXACT THE TRUE SOLUTION VALUES 



XMAX 

EPS 

MAX 

YP 

YC 

FP 

H 

ICON 



ERROR 

FCT 

EXACT 

RKUTT A_ 



THE ABSOLUTE DEVIATION 
'SOLUTION FROM THE TRUE 
SUBROUTINE TO 



FUNCTION 

'FUNCTION 

FUNCTION 

SOLUTION 



OF THE NUMERICAL 
SOLUTION 
EVALUATE ODE 



SUBROUTINE 
VALUES 
THE SUBROUTINE USED TO 
NEEDED STARTING VALUES 

» 4. »l» 4. 4* 4' 4# 4/ 4» 4- 4» 4. 4<- sV 4* 
— * — — .p. -r n~ n' n- nr -v> 



TO COMPUTE TRUE 



GENERATE 

- 4# 4* 4> 4/ 4« 4^ 4» nX- 4. 4» 4# 4. 4. . 

''r<rv¥^vv'i'T¥¥'rv^' ; r' 



INITIALIZE INPUT CONSTANTS 

XMAX=1 0.0 
MAX=2 

EPS=0. 000001 
K=3 

READ STEP SIZE AND CONTROL PRINTING 

READ! 5 , 100) H,ICON 

SPECIFY INITIAL CONDITION 

X0=0.0 
YO=l .0 
YR=Y0 

TEST FOR END OF INPUT DATA 

IF ( I CON— 0 ) 45,45,7 

PRINT INPUT DATA SET 

WR IT E ( 6 , 1000 ) H, ICON 
WRITE (6 ,1001 ) XO , YO 

COMPUTE NEXT SOLUTION POINT 
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if if if if if Me if Mr -5$ if Mr if -if if Mr if Mr M if- if if if if if if if if if if 



non non non a no non non non 



NC= ( XMAX+H/2. ) /H 
X1=X0+H 
X2=X0+2. *H 
X3=X0+3.*H 
X=X0+4.*H 

CALL RKUTTA{XO,YR, Y1,Y2,Y3,K,H) 

F1=FCT ( XI ,Y1 ) 

F2=FCT(X2,Y2) 

F3=FCT (X3,Y3) 

00 40 N=4 » NC 
M= 1 

YP=Y0+4.*H/3.*(2.*F1-F2+2.*F3) 

FP=FCT ( X » YP) 

11 YC= Y2+H/3 . *(F2+4.#F3+FP) 

DELY=ABS ( YC-YP ) 

TEST FOR CONVERGENCE 

I F { DELY-EPS ) 30,30,15 

TEST IF MAXIMUM NUMBER OF ITERATIONS IS SATISFIED 

15 IF ( M-MAX) 20,20,30 
20 FP=FCT ( X, YC ) 

M=M+1 
GO TO 11 

COMPUTE TRUE SOLUTION VALUE 
30 YEXACT=EXACT ( X) 

COMPUTE ERROR IN CALCULATION 
ERR0R=YEXACT-YC 

TEST IF DESIRED POINT OF PRINTING IS REACHED 

I F { (N/ICON*ICQN).EQ.N). WRITE (6, 1002) X , YP , YC , YEXACT , 
1ERR0R 

UPDATE REQUIRED PARAMETER VALUES 

F1=F2 
F2 = F 3 
F3=FP 
Y0=Y 1 
Y1=Y 2 
Y2=Y3 
Y3=YC 
X=X+H 

40 CONTINUE 

LOOP BACK TO READ NEXT SET OF INPUT DATA 

GO TO 5 
45 STOP 

100 FORMAT {F10.7 , 1 10) 

1000 FORMAT (////41X,'*INPUT DATA SET US ED* ' // 4 3X , • H =', 

1F10.7/43X, * ICON = * ,110) 

1001 FORMAT (//29X, '^RESULTS OF MILNE PREDICTOR-CORRECTOR 
1METH0D**'//10X, • X »,5X,« YP *,5X,' YC ', 
15X, • YEXACT ' , 5X , ' ERROR • //12X, F 1 0. 7 ,44X , 
1F10.7) 

1002 FORMAT ( '0',9X,F15.7,2X,F15.7,2X,F15.7,2X,F15.7,2X, 
1F15.7) 

END 
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SUBROUT INE RKUTT A ( XX , YY , Y1 , Y2 , Y3 , KK , HH ) 

NN=1 

17 CK1=HH*FCT( XX, YY) 

CK2=HH*FCT (XX+HH/2.0, YY+CK 1/2.0) 

CK3=HH*FCT ( XX+HH /2 • 0 , YY+CK2/2 . 0 ) 
CK4=HH*FCT(XX+HH, YY+CK3) 

GO TO (101,102,103) ,NN 

101 Yl=YY+(CKl+2.0*CK2+2.0*CK3+CK4)/6. 0 
YY = Y 1 

77 XX=XX+HH 

YEXACT=EXACT(XX) 

ERR0R= Y EXACT- YY 

WRITE(6 ,1004 ) XX ,YY,YEXACT , ERROR 
NN=NN+ 1 

IF ( iMN-KK > 17,17,99 

102 Y2 = YY+ lCKl+2.0*CK2+2.0*CK3+CK4)/6.0 
YY=Y2 

GO TO 77 

103 Y3=YY+ (CK1+2.0*CK2+2.0*CK3+CK4 )/6.0 
YY=Y3 

GO TO 77. 

1004 FORMAT ( »0',9X,F15.7,2X,F15.7,16X,F15.7,2X,F15.7) 
99 RETURN 
END 



FUNCTION FCT ( XN , YN) 

FCT=— Y N 

RETURN 

END 



FUNCTION EXACT (XK ) 

EXACT=EXP(-XK) 

RETURN 

END 
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* INPUT DATA SET USED* 

H = 0.5000000 
ICON = 2 

**RESULTS OF MILNE PREDICTOR-CORRECTOR METHOD** 



X 


YP 


YC 


YEXACT 


ERROR 


0.0000000 






1.0000000 




0. 5000000 


0.6067709 




0.6065306 


-0.0002403 


1.0000000 


0.3681710 




0.367894 


-0.0002916 


1.5000000 


0.2233955 




0.2231301 


-0.0002654 


2. 0000000 


0.1385594 


0.1353100 


0.1353353 


0.0000253 


3.0000000 


0.0509259 


0.0496315 


0.0497871 


0.0001555 


4. 0000000 


0.0183159 


0.0181091 


0.0183156 


0.0002066 


5.0000000 


0.0061856 


0.0064749 


0.0067379 


0.0002630 


6.0000000 


0.0015049 


0. 0021262 


0.0024788 


0.0003525 


7. 0000000 


-0.0005340 


0.0004233 


0.0009119 


0.0004886 


8.0000000 


-0.0017356 


-0. 0003513 


0.0003355 


0.0006868 


9. 0000000 


-0. 0028179 


-0.0008470 


0.0001234 


0.0009704 


10.0000000 


-0.0041233 


-0.0013280 


0.0000454 


0.0013734 



* I NPU T DATA SET USED* 

H = 1.0000000 
ICON = 1 

**RESULT S OF MILNE PREDICTOR-CORRECTOR METHOD** 



X 

0.0000000 

1.0000000 

2.0000000 

3.0000000 

4. 0000000 

5.0000000 

6. 0000000 

7.0000000 

8.0000000 

9.0000000 

10.0000000 



YP 

0.3750000 

0.1406250 

0.0527344 

0.0468752 

0.0147570 

0.0102883 

0.0014462 

-0.0005506 

-0.0039288 

-0.0064989 



YC 



0.0164931 
0.0051922 
0.0002441 
0.0005433 
-0. 0009222 
0.0010040 
-0.0017257 



YEX ACT 

1.0000000 
0.3678794 
0.1353353 
0.0497871 
0.0183156 
0.0067379 
0.0024788 
0.0009119 
0.0003355 
0.0001234 
0.0000454 



ERROR 

-0.0071206 

-0.0052897 

-0.0029473 

0.0018226 

0.0015457 

0.0022346 

0.0003686 

0.0012576 

-0.0008806 

0.0017711 
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PROGRAM THREE 



NYSTROM PREDICTOR-CORRECTOR METHOD. ODE: Y'=-Y 
WIYH Y( 0 ) = 1 . TO SOLVE ANOTHER ODE SIMPLY CHANGE 
INITIAL CONDITIONS AND FUNCTION SUBROUTINES. 



C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



c 

c 

c 

5 

c 

c 

c 



c 

c 

c 

c 

c 

c 

7 

c 

c 

c 



%*<■ %>. yu v. VL ^ ou ^ vu j. <a> ^ J* «ju J. •*. ^ <ju <ju «ju j. yu ou 

^ ^ ^ ^ ^ r^s #j% ^ nj" rf* -'p #p *p *p -p ^ «p 4|« «p «fp #p 



* FORTRAN PROGRAM FOR NYSTROM PREDICTOR-CORRECTOR 

* METHOD IN THE NUMERICAL SOLUTION OF THE ORDINARY 

* DIFFERENTIAL EQUATION IN THE FORM OF DY/DX=F(X,Y) 

* WITH THE INITIAL CONDITION Y(XO)=YO 

* THE PARAMETERS TO THE PROGRAM HAVE THE FOLLOWING 

* MEANINGS... 

* XMAX THE TERMINAL BOUNDARY CONDITION 

* EPS_ THE CONVERGENCE TEST CONSTANT 

* MAX THE MAXIMUM NUMBER OF ITERATIONS 

* YP THE PREDICTED VALUE OF THE NEXT 

* SOLUTION POINT 

* YC_, THE CORRECTED VALUE OF THE NEXT 

* SOLUTION POINT 

* FP THE FUNCTION EVALUATION 

* H THE STEP SIZE TO BE USED TO ADVANCE 

* “THE POINT OF INTEGRATION 

* ICON THE DESIRED POINT OF PRINTING THE 

* COMPUTED SOLUTION FOR EVERY FIXED 

* NUMBER OF INTEGRATION POINTS 

* Y EXACT THE TRUE SOLUTION VALUES 

* ERROR THE ABSOLUTE DEVIATION OF THE NUMERICAL 

* SOLUTION FRGM THE TRUE SOLUTION 

* FCT FUNCTION SUBROUTINE TO EVALUATE ODE 

* FUNCTION. 

* EXACT FUNCTION SUBROUTINE TO COMPUTE TRUE 

* SOLUTION VALUES 

* RKUTTA THE SUBROUTINE USED TO GENERATE 

* NEEDED STARTING VALUES 



INITIALIZE INPUT CONSTANTS 

XMAX= 10 . 0 
MAX=2 

EPS=0. 000001 

READ STEP SIZE AND CONTROL PRINTING 
READ 15,100) H , I C ON 
SPECIFY INITIAL CONDITION 



X0=0.0 
Y0=1 .0 
YR = YO 

TEST FOR END OF INPUT DATA 

I F ( I CON— 0 1 45,45,7 

PRINT INPUT DATA SET 

WRIT E ( 6 , 1000 ) H , ICON 
WRITE! 6,1001) XO , YO 



POINT 



COMPUTE NEXT SOLUTION 
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it # if -tf it it it it * it it it it it it * it it it it it it- it if it- if -K- if it it 



ooo non ooo non ono noo ooo 



NC= ( XMAX+H/2. ) /H 
X 1=X0+H 
X=X0+2.*H 

CALL RKUTTA(X0,YR,Y1,H) 

F1=FCT (XI, Yl) 

DO 40 N=2,NC 
M= 1 

YP=Y0+2.*H^F1 
FP=FCT (X,YP) 

11 YC=Y1+H/2.*(F1+FP) 

DELY=ABS ( YC-YP ) 

TEST FOR CONVERGENCE 

IF (DELY-EPS) 30,30*15 

TEST IF MAXIMUM NUMBER OF ITERATIONS IS SATISFIED 

15 IF { M-MAX ) 20,20,30 

20 FP=FCT ( X , YC ) 

M=M+ 1 
GO TO 11 

COMPUTE TRUE SOLUTION VALUE 
30 Y EX ACT = EXACT ( X ) 

COMPUTE ERROR IN CALCULATION 
ERR0R=YEX ACT-Y C 

TEST IF DESIRED POINT OF PRINTING IS REACHED 

IF ( ( N/ICON* ICON) .EQ.N) WR I T E (6 , 1 002 ) X , YP, YC , YEXACT , 
1ERR0R 

UPDATE REQUIRED PARAMETER VALUES 

Y0=Y 1 
Y1=YC 
F1=FP 
X = x + H 

40 CONTINUE 

LOOP BACK TO READ NEXT SET OF INPUT DATA 

GO TO 5 
45 STOP 

100 FORMAT (F10-7, 110) 

1000 FORMAT (////41X,**INPUT DATA SET US ED* • //43X , • H =', 

1F10.7/43X. • ICON = * ,110) 

1001 FORMAT (// 29 X,'* RESULTS OF NYSTROM PREDICTOR-CORRECTOR 
1METH0D *'//10X,' X • ,5X,‘ YP ',5X,« YC ', 

15X , ' YEXACT * , 5X , ' ERROR • //12X, F10. 7 ,44X , 

1F10.7) 

1002 FORMAT { * O' , 9X , F 1 5. 7 , 2X*F 15 . 7 ,2 X , FI 5 . 7 , 2X , F15 . 7 , 2X , 
1F15.7) 

END 



SUBROUTINE RKUTT A ( XX , YY , Yl , HH ) 

CK1=HH*FCT(XX,YY) 

CK2=HH*FCT( XX+HH/2.0,Y Y+CK 1/2.0) 

CK3=HH*FCT ( XX+HH/2 . 0 , YY+CK2 /2 i 0 ) 

CK4=HH*FCT (XX+HH ,YY+CK3) 

Yl=YY+(CKl+2.0*CK2+2. 0*CK3+CK4) /6. 0 

XX=XX+HH 

YEXACT=EXACT (XX) 

ERROR=Y EXACT- Y 1 

WRITE (6, 10041 XX, Yl, YEXACT , ERROR 
1004 FORMAT ( , 0',9X,F15.7,2X,F15.7,16X,F15.7,2X,F15.7) 
RETURN 
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END 



FUNCTION FCT(XN,YN) 

FCT=-YN 

RETURN 

END 



FUNCTION EXACT { X K ) 
EXACT=EXP (-XK) 
RETURN 
END 
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* INPUT DATA SET USED* 

H = 0.5000000 
ICON = 2 

^RESULTS OF NYSTROM PREDICTOR-CORRECTOR METHOD * 



X 


YP 


YC 


YEXACT 


ERROR 


0.0000000 






1.0000000 




0. 5000000 


0.6067709 




0.6065306 


-0.0002403 


1.0000000 


0.3932291 


0 .3636068 


0.3678794 


0.0042726 


2.0000000 


0.1444499 


0. 1298206 


0.1353353 


0.0055147 


3.0000000 


0.0516075 


0.0463004 


0 .0497871 


0.0034867 


4.0000000 


0.0184066 


0.0165119 


0.0183156 


0.0018037 


5.0000000 


0.0065643 


0.0058886 


0.0067379 


0.0008494 


6.0000000 


0.0023410 


0.0021000 


0.0024788 


0.0003788 


7.0000000 


0.0008349 


0. 0007489 


0.0009119 


0.0001630 


8.0000000 


0.0002977 


0.0002671 


0 .0003355 


0.0000684 


9.0000000 


0.0001062 


0.0000952 


0.0001234 


0.0000282 


10.0000000 


0.0000379 


0.0000340 


0.0000454 


0.0000114 



* INPUT DATA SET USED* 

H = 1.0000000 
ICON = 1 

*RESULT S OF NYSTROM PREDICTOR-CORRECTOR METHOD * 

X YP YC YEXACT ERROR 



0.0000000 






1.0000000 




1.0000000 


0.3750000 




0.3678794 


-0.0071206 


2.0000000 


0.2500000 


0. 1093750 


0.1353353 


0.0259603 


3.0000000 


0. 0625000 


0.0156250 


0.0497871 


0.0341621 


4.0000000 


0.0468750 


-0.0053594 


0.0183156 


0.0241750 


5.0000000 


-0.0078125 


-'0. 0078125 


0.0067379 


0.0145504 


6.0000000 


0.0097656 


-0.0041504 


0.0024788 


0.0066291 


7.0000000 


-0.0087891 


-0.0021973 


0.0009119 


0. 0031091 


8.0000000 


0. 0046387 


-0.0005798 


0.0003355 


0.0009153 


9.0000000 


-0.0045166 


-0 .0003052 


0.0001234 


0.0004286 


10.0000000 


0.0028381 


0.0000572 


0.0000454 


-0.0000118 
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PROGRAM FOUR 



HERMITE PREDICTOR-CORRECTOR METHOD. ODE: Y»=-Y 
WIYH Y( 0 ) =1 . TO SOLVE ANOTHER ODE SIMPLY CHANGE 
INITIAL CONDITIONS AND FUNCTIGN SUBROUTINES. 



XMAX 

EPS 

MAX 

YP_ 



YC 

FP 

H 

ICON 



YE X ACT 
ERROR 



THE 




4^ Vf ^ Vf ^ ^ %v v? ^ ^ ^ ^ ^ ^ v# ^ ^ ^ ^ ^ ^ u# y** %v ^ vL %u 

^ ^ ^ ^ ^p >p ^p **p <rj% #p #p #p -'p -*p ^p *p pp fp ^p /p -p *p «p ^p ^p» #p #w ^p 

* FORTRAN PROGRAM FOR HERMITE PREDICTOR-CORRECTOR 

* METHOD IN THE NUMERICAL SOLUTION OF THE ORDINARY 

* DIFFERENTIAL EQUATION IN THE FORM OF DY/DX=F(X,Y) 

* WITH THE INITIAL CONDITION Y(X0)=Y0 

* THE PARAMETERS TO THE PROGRAM HAVE THE FOLLOWING 
MEANINGS. 

TERMINAL BOUNDARY CONDITION 
CONVERGENCE TEST CONSTANT 
MAXIMUM NUMBER OF ITERATIONS 
PREDICTED VALUE OF THE NEXT 
SOLUTION POINT 

THE CORRECTED VALUE OF THE NEXT 

SOLUTION POINT 

THE FUNCTION EVALUATION 

THE STFP SIZE TO BE USED TO ADVANCE 

THE POINT OF INTEGRATION 
THE DESIRED POINT OF PRINTING THE 
COMPUTED SOLUTION FOR EVERY FIXED 
NUMBER OF INTEGRATION POINTS 
THE TRUE SOLUTION VALUES 

THE ABSOLUTE DEVIATION OF THE NUMERICAL 
SOLUTION FROM THF TRUE 
SUBROUTINE TO 



* FCT FUNCTION 

* FUNCTION 

* EXACT FUNCTION 

* SOLUTION 



SOLUTION 

EVALUATE 



ODE 



SUBROUTINE 
VALUES 

RKUTTA THE SUBROUTINE USED TO 

NEEDED STARTING VALUES 



TO COMPUTE TRUE 



GENERATE 



Vf ^ V* V? V* V# Jf V? Vr V' V* U# JU V* V* V* Vf ^ VU V^ V* ^ ^ V* W V# ^ ^ V* V^ ^ % 

^P ^p #p #p ^ >p -p /p pp ^p <Y* ^ ^p ^p pp «^p #p ^p ^p *p #p 4p #p #p ^p ^ ^ /p -V- /p ^ r % >p -*p ^ >p ^ V* ■"** # 



5 



7 



INITIALIZE INPUT CONSTANTS 

XM AX=1 0.0 
MAX= 2 

EPS=0. 000001 

READ STEP SIZE AND CONTROL PRINTING 

RE AD (5*100) H* I CON 

SPECIFY INITIAL CONDITION 

X0=0. 0 
YO= 1 . 0 
YR=Y0 



TEST FOR END OF INPUT DATA 

I F ( I CON-O ) 45,45,7 

PRINT INPUT DATA SET 

WRITE( 6,1000) H , I C ON 
WRITE(6, 1001 ) XO , Y 0 



POINT 



COMPUTE NEXT SOLUTION 
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NC=(XMAX+H/2.)/H 
Xl = XO+H 
X=X0+2.*H 

CALL RKUTTA(X0,YR,Y1 ,H) 

FO=FCTl XO, YO) 

F1 = FCT { X 1 , Y1 ) 

DO 40 N=2,NC 
M=1 

YP=5.*Y0-4.*Y1+H*(4.*F1+2.*F0) 

FP=F CT ( X , Y P) 

11 YC=Y0+H/3.*(F0+4.*F1+FP) 

DELY=ABS(YC-YP) 

TEST FOR CONVERGENCE 

IF ( DELY-EPS) 30,30,15 

TEST IF MAXIMUM NUMBER OF ITERATIONS IS SATISFIED 

15 IF(M-MAX) 20,20,30 
20 FP= FCT ( X , YC ) 

M=M+1 
GO TO 11 

COMPUTE TRUE SOLUTION VALUE 
30 Y EXACT = EXACT ( X ) 

COMPUTE ERROR IN CALCULATION 
ERROR=YEXACT-YC 

TEST IF DESIRED POINT OF PRINTING IS REACHED 

IF ( ( N/ ICON* ICON ) . EQ.N ) WRI T E { 6 , 1002 ) X, YP, YC , YEXACT, 
1ERR0R 

UPDATE REQUIRED PARAMETER VALUES 

Y0=Y 1 
Y 1 =Y C 
F0=F 1 
F 1= F P 
X=X + H 

40 CONTINUE 

LOOP BACK TO READ NEXT SET OF INPUT DATA 

GO TO 5 
45 STOP 

100 FORMAT (F10.7 , 1 10) 

1000 FORMAT 1////41X, ' *1 NPUT DATA SET US ED* ' //43X , ' H =', 
1 F 1 C . 7 / 43 X , ' ICON =* , 110) 

1001 FORMAT (//29X,'* RESULTS OF HERMITE PREDICTOR-CORRECTOR 
1METH0D **//10X,' X ',5X,» YP ',5X,' YC ', 

1 5X , ' YEXACT ' , 5X , ' ERROR ' // 12 X, F 1 0 . 7 , 44X , 

1F10. 7) 

1002 FORMAT { *0',9X,F15.7,2X,F15.7,2X,F15.7,2X,F15.7,2X, 
1F15.7) 

END 



SUBROUTINE RKUTT A( XX , YY , Y1 , HH ) 
CK1=HH*FCT(XX, YY) 

CK2=HH*FCT(XX+HH/2.0, YY+CK 1/2.0) 

CK3=HH*FCT (XX+HH/2.0 ,YY+CK2/2.0) 

CK4=HH*FCT( XX+HH , YY+CK3) 

Yl=YY+(CKl+2.0*CK2+2.0*CK3+CK4)/6.0 

XX=XX+HH 

YEXACT=EXACT(XX) 

ERR0R=YEXACT-Y1 

WR I T E ( 6 , 1 004 ) XX, Yl, YEXACT, ERROR 
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1004 FORMAT ( «0',9X,F15.7.2X»F15.7»16X,F15.7,2X,F15.7) 
RETURN 
END 



FUNCTION FCT( XN* YN) 

FCT=-YN 

RETURN 

END 



FUNCTION EXACT { X K ) 
EXACT=EXP (— XK) 
RETURN 
END 
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* INPUT DATA SET USED* 

H = 0.5000000 
ICON = 2 

^RESULTS OF HERMITE PREDICTOR-CORRECTOR METHOD * 



X 


YP 


YC 


YEXACT 


ERROR 


0.0000000 






1.0000000 




0.5000000 


0.6067709 




0.6065306 


-0.0002403 


1.0000000 


0.3593760 


0.3675976 


0.3678794 


0. 0002818 


2. 0000000 


0. 1296880 


0.1349390 


0.1353353 


0.0003963 


3.0000000 


0.0439692 


0.0491728 


0.0497871 


0.0006143 


4.0000000 


0.0106042 


0. 0173813 


0.0183156 


0.0009343 


5.0000000 


-0.0043826 


0.0053371 


0.0067379 


0.0014008 


6.0000000 


-0.0139239 


0.0003908 


0.0024788 


0.0020879 


7. 0000000 


-0.0234190 


-0.0021938 


0.00091 19 


0.0031057 


8.0000000 


-0.0358057 


-0.0042811 


0.0003355 


0.0046165 


9.0000000 


-0.0535792 


-0.0067376 


0.0001234 


0.0068610 


10.0000000 


-0.0797584 


-0.0101509 


0.0000454 


0.0101963 



*INPUT DATA SET USED* 

H = 1.0000000 
ICON = 1 

*RESULTS OF HERMITE PREDICTOR-CORRECTOR METHOD * 



X 


YP 


YC 


YEXACT 


ERROR 


0. 0000000 






1.0000000 




1.0000000 


0.3750000 




0.3678794 


-0.0071206 


2.0000000 


0.0000000 


0.1296299 


0.1353353 


0.0057054 


3.0000000 


0.1620350 


0.0732166 


0.0497871 


-0.0234295 


4.0000000 


-0.2105594 


-0.0092715 


0.0183156 


0.0275871 


5.0000000 


0.3834665 


0. 0599074 


0 .0067379 


-0.0531695 


6.0000000 


0.6344536 


-0.0839149 


0.0024788 


0.0863936 


7.0000000 


1.0731880 


0. 1479157 


0.0009119 


-0.1470038 


8.0000000 


-1.8065000 


-0.2466852 


0 .0003355 


0.2470207 


9.0000000 


3.0441850 


0.4165850 


0.0001234 


-0.4164615 


10.0000300 


-5.1285590 


-0.7014816 


0 .0000454 


0.7015270 
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PROGRAM FIVE 



HAMMING PREDICTOR-CORRECTOR METHOD. ODE: Y'=-Y 
WIYH Y( 0 ) = 1 . TO SOLVE ANOTHER ODE SIMPLY CHANGE 
INITIAL CONDITIONS AND FUNCTION SUBROUTINES. 



C 

C 

C 

C 

C 

C 

c 



V- O- J# ^ Ur ^ ^ *4- >0 U/ ^ V# Uo X -JU mJr %t 4# ^ 4/ ^ ^ J# V# Jf X ^IU 

••j** 'Y' n* y> *y* i 4 v »{% «y* ^ ^ ^ *y» ■'"p ^ «y* 'P ^p ^p- #p «y» «y* #;• ip #p «p «y* «p #p pp ip ^ 

FORTRAN PROGRAM FOR HAMMING PREDICTOR-CORRECTOR 



C 

c 

c 

c 

c 



c 

c 

c 

c 

c' 

c 

c 



c 

c 

c 

c 

c 

c 



* 

* 

* 

❖ 



METHOD IN THE NUMERICAL SOLUTION OF THE ORDINARY 
DIFFERENTIAL EQUATION IN THE FORM OF DY/DX=F(X,Y) 
WITH THE INITIAL CONDITION Y(XO)=YO 
THE PARAMETERS TO THE PROGRAM HAVE THE FOLLOWING 
MEANINGS... 

THE TERMINAL BOUNDARY CONDITION 
E CONVERGENCE TEST CONSTANT 
THE MAXIMUM NUMBER OF ITERATIONS 
EDICTED SOLUTION VALUE 
MODIFIED PREDICTED SOLUTION VALUE 
CORRECTED SOLUTION VALUE 
MODIFIED CORRECTED SOLUTION VALUE 
THE FUNCTION EVALUATION 
THE STEP SIZE TO BE USED 
THE POINT OF INTEGRATION 

E DESIRED POINT OF PRINTING THE 
COMPUTED SOLUTION FOR EVERY FIXED 
NUMBER OF INTEGRATION POINTS 
E TRUE SOLUTION VALUES 



c 


* 


XMAX 


c 


"T* 


EPS 


c 


£ 


MAX 


c 


* 


YPM 


c 


J. 

-r* 


YP i 


c 


* 


YCffi 


c 


* 


YC 


c 


* 


FP 


c 




h ; 


c 


* 




c 


* 


ICON 


c 


* 


1 


c 


* 


1 


c 




YEX ACT 


c 




ERROR 


c 


* 




c 


* 


FCT 


c 


❖ 




c 


* 


EXACT 


c 


❖ 




c 

r 


* 

WLr 


RKUTTA 

1 



TO ADVANCE 



ABSOLUTE DEVIATION 
SOLUTION FROM THE TRUE 
SUBROUTINE TO 



FUNCTION 



OF THE NUMERICAL 
SOLUTION 
EVALUATE ODE 



TO COMPUTE TRUE 



SUBROUTI NE 

SOLUTION VALUES 

THE SUBROUTINE USED TO GENERATE 
NEEDED STARTING VALUES 

%U vV *V 'V ***** V 4 * 4# %l*> ^ «JU JU Vf v 1 ^ JU J# Oy v** J# JU «JU %'# J/ %V >Lr JU JU JU wu wu 

ip «p #p ip ^ ip ip ip ip ^P ip ^p ip ^p ip i<fc ip ip ^p #P ^p **P *Y^ "Ip ^P #p ip ip «Y* ^p ^ rp ip ijk *p i|S ^p ip #p ip Kp rp ^p 



INITIALIZE INPUT CONSTANTS 

XMAX= 10.0 
MAX=2 

EPS=0. 000001 
K=3 

READ STEP SIZE AND CONTROL PRINTING 

READ ( 5 1 100 ) H, I CON 

SPECIFY INITIAL CONDITION 

ET-0.0 
X0=0. 0 
Y0= 1 .0 
YR=Y0 

TEST FOR END OF INPUT DATA 

I F ( I CON— 0 ) 45,45,7 

PRINT INPUT DATA SET 

WRITEI6, 1000) H, ICON 
WRITE(6,1001 ) XO , Y 0 
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ono ooo ooo non ooo ooo non on 



COMPUTE NEXT SOLUTION POINT 



NC=( XMAX+H/2. ) /H 

X1=X0+H 

X2=X0+2.*H 

X3=X0+3.*H 

X=X0+4.*H 

CALL RKUTTA(XO,YR, Y1,Y2,Y3,K,H) 

F1=FCT(X1,Y1) 

F2=FCT ( X2, Y2) 

F3=FCT (X3,Y3) 

DO 40 N=4 1 NC 
M= 1 

YPM=Y0+4.*H/3.*(2.*F1-F2+2.*F3) 

YP=YPM-112./121.*ET 
FP=FCT ( X , YP ) 

11 YCM=1./8.*(9.*Y3-Y1 )+3.*H/8.*( FP+2«*F3-F2) 

ET= YPM— YCM 

YC=YCM+9./121.*ET 

DELY=ABS(YC--YP) 

TEST FOR' CONVERGENCE 

I F { DEL Y-EP S) 30,30,15 

TEST IF MAXIMUM NUMBER OF ITERATIONS IS SATISFIED 

15 IF ( M— MAX ) 20,20,30 

20 FP=FCT (X,YC) 

M=M+ 1 
GO TO 11 

COMPUTE TRUE SOLUTION VALUE 
30 YEXACT=EXACT { X ) 

COMPUTE ERROR IN CALCULATION 
ERROR=YEXACT-YC 

TEST IF DESIRED POINT OF PRINTING IS REACHED 

I F ( ( N/ ICON* ICON). E Q. N) WRI T E (6 , 1 00 2 ) X ,YP , YC , YEXACT , 
1 ERROR 

UPDATE REQUIRED PARAMETER VALUES 



F3=FP 
Y0=Y 1 
Y1=Y2 
Y2=Y3 
Y3=Y C 
X — X + H 

40 CONTINUE 

LOOP BACK TO READ NEXT SET OF INPUT DATA 

GO TO 5 
45 STOP 

100 FORMAT (F10.7 ,1 10) 

1000 FORMAT (////41X, ' * I N P UT DATA SET US ED* • //43X , ' H =', 
1F10.7/43X, • ICON =',1101 

1001 FORMAT! //29X, ' *RESULTS OF HAMMING PREDICTOR-CORRECTOR 
1METH0D *‘//10X,' X ',5X,' YP ',5X,' YC ', 

1 5X , ' YEXACT ' , 5X , ' ERROR • / / 12X , F 1 0 . 7 , 44X , 

1F10. 7) 

1002 FORMAT ( • 0 ' ,9 X , F 1 5 . 7 , 2X , F 15 . 7 , 2X, FI 5 . 7 , 2X , F 1 5 . 7 , 2X , 
1F15.7) 

END 
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SUBROUTINE RKUTT A ( XX , Y Y , Y1 , V2 , V3 , KK , HH ) 

NN=1 

CK1=HH*FCT (XX, YY ) 

CK2=HH* t FCT (XX+HH/2 .0,YY+CK 1/2.0) 
CK3=HH*FCT(XX+HH/2.0,YY+CK2/2.0) 

CK4=HH*FCT (XX+HH, YY+CK3 ) 

GO TO ( 101 ,102, 103 ) ,NN 

101 Yl=YY+(CKl+2.0*CK2+2.0*CK3+CK4)/6. 0 
YY=Y 1 

77 XX=XX+HH 

YEXACT=EXACT (XX) 

ERROR= YEXACT-YY 

WRITE! 6,1004) XX , YY , YE XACT , ERROR 
NN=NN+ 1 

IF(NN-KK) 17,17,99 

102 Y2=YY+(CKl+2.0*CK2+2.0*CK3+CK4)/6. 0 
YY=Y2 

GO TO 77 

103 Y3=YY+(CKl+2.0*CK2+2.0*CK3+CK4)/6. 0 
YY=Y3 

GO TO 77 

1004 FORMAT! ‘0',9X,F15.7,2X,F15.7,16X,F15.7,2X,F15.7) 
99 RETURN 
END 



FUNCTION FCT ( XN , YN ) 

FCT=-Y N 

RETURN 

END 



FUNCTION EXACT(XK) 
EXACT=EXP (— XK) 
RETURN 
END 
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* I NPUT DATA SET USED* 

H = 0.5000000 
ICON = 2 

^RESULTS OF HAMMING PREDICTOR-CORRECTOR METHOD * 



X 


YP 


YC 


YEXACT 


ERROR 


0.0000000 






1.0000000 




0.5000000 


0.6067709 




0.6065306 


-0.0002403 


1.0000000 


0.3681710 




0.3678794 


-0.0002916 


1.5000000 


0.2233955 




0.2231301 


-0.0002654 


2.0000000 


0.1385594 


0.1355409 


0.1353353 


-0. 0002056 


3. 0000000 


0. 0494538 


0. 0499289 


0.0497871 


-0.0001418 


4.0000000 


0.0184383 


0.0183963 


0.0183156 


-0.0000806 


5.0000000 


0.0068499 


0.0067794 


0.0067379 


-0.0000415 


6.0000000 


0.0025701 


0.0024990 


0.0024788 


-0.0000202 


7.0000000 


0.0009748 


0. 0009215 


0.0009119 


-0.0000096 


8.0000000 


0.0003766 


0.0003401 


0.0003355 


-0.0000046 


9.0000000 


0.0001496 


0.0001256 


0.0001234 


-0.0000022 


10.0000000 


0.0000619 


0.0000465 


0.0000454 


-0.0000011 





*INPUT 


DATA SET USED* 






H 


= 1.0000000 








ICON 


= 1 






^RESULTS OF HAMMING 


PREDICTOR- 


CORRECTOR 1 


METHOD * 


X 


YP 


YC 


YEXACT 


ERROR 


0.0000000 






1.0000000 




1.0000000 


0.3750000 




0 .3678794 


-0.0071206 


2.0000000 


0. 1406250 




0.1353353 


-0. 0052897 


3.0000000 


0. 0527344 




0.0497871 


-0.0029473 


4.0000000 


0.0468752 


0.0190868 


0.0183156 


-0.0007712 


5.0000000 


-0.0199183 


0.0056581 


0.0067379 


0.0010798 


6. 0000000 


0.0245465 


0.0057380 


0.0024788 


-0.0032592 


7.0000000 


-0.0516121 - 


0. 0008928 


0.0009119 


0. 0018047 


8.0000000 


0. 0793315 


0.0053606 


0.0003355 


-0.0050251 


9.0000000 


-0.1184798 - 


0.0069360 


0.0001234 


0.0070594 


10.0000000 


0.1838089 


0. 0110968 


0.0000454 


-0.0110514 
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PROGRAM SIX 



SECOND ORDER ADAMS PREDICTOR-CORRECTOR METHOD. ODE: 
Y • =— Y WITH Y { 0 ) = 1 . TO SOLVE ANOTHER ODE CHANGE 
INITIAL CONDITIONS AND FUNCTION SUBROUTINES. 



***** 



O. ^ V ' «JU 0<» X J# Jf J** >■> <JU sAf V" J. J. V* 'A- «*» <JU OU J. 

V 'i' V v v v 'i' r <i % VT'p r 'i' t t r r v *1 ■ r ¥ t V ' • V t r t t *1' J i- *i ■ t r « r ¥ 'r' *r 



* FORTRAN PROGRAM FOR SECOND ORDER ADAMS PREDICTOR— 

* CORRECTOR METHOD IN THE NUMERICAL SOLUTION OF THE 

* ORDINARY DIFFERENTIAL EQUATION IN THE FORM OF 

* DY/DX=F(X,Y) WITH THE INITIAL CONDITION Y(X0)=Y0 

* THE PARAMETERS TO THE PROGRAM HAVE THE FOLLOWING 

* MEANINGS... 

* XMAX THE TERMINAL BOUNDARY CONDITION 

* EP S__ THE CONVERGENCE TEST CONSTANT 

* MAX THE MAXIMUM NUMBER OF ITERATIONS 

* YP THE PREDICTED VALUE OF THE NEXT 

* SOLUTION POINT 

* YC THE CORRECTED VALUE OF THE NEXT 

* SOLUTION POINT 

* FP THE FUNCTION EVALUATION 

* H THE STEP SIZE TO BE USED TO ADVANCE 

* THE POINT OF INTEGRATION 

* ICON. THE DESIRED POINT OF PRINTING THE 

* COMPUTED SOLUTION FOR EVERY FIXED 

* NUMBER OF INTEGRATION POINTS 

* YEXACT THE TRUE SOLUTION VALUES 

* ERROR THE ABSOLUTE DEVIATION OF THE NUMERICAL 

* SOLUTION FROM THE TRUE SOLUTION 

* FCT FUNCTION SUBROUTINE TO EVALUATE ODE 

* FUNCTION, 

* EXACT FUNCTIGN SUBROUTINE TO COMPUTE TRUE 

* SOLUTION VALUES 

* RKUTTA THE SUBROUTINE USED TO GENERATE 

* NEEDED STARTING VALUES 



5 



7 



INITIALIZE INPUT CONSTANTS 

XMAX=1 0 . 0 
MAX=2 

EPS=0. OOOOOl 

READ STEP SIZE AND CONTROL PRINTING 

READ (5 » 100 ) H» ICON 

SPECIFY INITIAL CONDITION 

X0=0. 0 
Y0= 1 . 0 
YR=Y0 



TEST FOR END OF INPUT DATA 

I F { I CON-O ) 45 1 45 1 7 

PRINT INPUT DATA SET 

WRITE ( 6, 1000) H » I C ON 
WRITEI6, 1001) XO i YO 



POINT 



COMPUTE NEXT SOLUTION 
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to to to to to to to to to to to to to to to to to to to to to to -if to to to to to to -:j 



non non onn non non oon ooo 



NC= ( XMAX+H/2 • ) / H 
FO=FCT(XO,YO) 

X1=X0+H 
X=X0+2 .*H 

CALL RKUTTA{ XO,YR, Y1 ,H) 

F1=FCT ( X 1 , Y1 ) 

DO 40 N=2 » NC 
M= 1 

YP=Y1+H/2.*(3.*F1-F0) 

FP=FCT ( X »YP) 

11 YC=Y1+H/2.*(FP+F1) 

DELY=ABS ( YC-YP ) 

TEST FOR CONVERGENCE 

IF(DELY-EPS) 30,30,15 

TEST IF MAXIMUM NUMBER OF ITERATIONS IS SATISFIED 

15 I F ( M— MAX) 20,20,30 

20 FP= FCT ( X , YC ) 

M=M+1 
GO TO 11 

COMPUTE TRUE SOLUTION VALUE 
30 YEXACT=EXACT(X) 

COMPUTE ERROR IN CALCULATION 
ERROR=YEXACT-YC 

TEST IF DESIRED POINT CF PRINTING IS REACHED 

I F ( ( N/ I CON* ICON ) . EQ • N ) WRITE (6, 1002) X, YP , YC , YEXACT , 
1ERROR 

UPDATE REQUIRED PARAMETER VALUES 



40 



Y0=Y 1 
Y1=Y C 
F 1 =F P 
X=X+H 
CONTINUE 

LOOP BACK TO READ NEXT SET OF INPUT DATA 



GO TO 5 
45 STOP 

100 FORMAT ( F 10 • 7 , 1 1 0 ) 

1000 FORMAT {////41X, • *INPUT DATA SET USED*' //43X, 'H 
1F10.7/43X, ' ICON = • , I 1 0 ) 

1001 FORMAT ( // 29X , ' **RESULTS OF ADAM2 PREDICTOR-CORRECTOR 



= ' , 



1METHOD** ' // 10X , ' 
YEXACT 



, 5X : 



X 



' » 5X , 



ERROR 



YP ' , 5X , ' YC 
' //12X,F10.7 ,44X , 



15X, ' 

1 FI 0 . 7 ) 

1002 FORMAT ( • O' ,9X , FI 5 . 7 , 2X , F15 .7 ,2 X , FI 5 .7 , 2X , FI 5 . 7 , 2X , 
1F15.7) 

END 



SUBROUTINE RKUTTAt XX, YY, Y1 ,HH) 

CK1=HH*FCT (XX, YY ) 

CK2 = HH*FCT ( XX+HH/2 . 0 , YY+CK1 /2 * 0 ) 

CK3=HH*FCT( XX+HH/2. 0,YY+CK2/2. 0) 

CK4=HH*FCT ( XX+HH , YY+CK 3 ) 

Y1 =YY+ ( C K 1+2 . 0*C K2+2 . 0 *CK3 +CK4 ) / 6 . 0 
XX=XX+HH 

YEXACT=EXACT (XX) 

ERROR= YEXACT- Y1 

WRITE (6, 1004) XX, Yl, YEXACT, ERROR 
1004 FORMAT ('0',9X,F15.7,2X,F15.7,16X,F15.7,2X,F15.7) 
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RETURN 

END 



FUNCTION FCT (XN , YN ) 

FCT=- YN 

RETURN 

END 



FUNCTION EXACT ( XK) 
EXACT=EXP (-XK) 
RETURN 
END 



* IN PUT DATA SET USED* 

H = 0.5000000 
ICON = 2 

**RE SULTS OF A DA M2 PREDICTOR-CORRECTOR METHOD** 



X 


YP 


YC 


YEXACT 


ERROR 


0.0000000 






1.0000000 




0. 5000000 


0.6067709 




0.6065306 


-0.0002403 


1.0000000 


0.4016927 


0.3634746 


0.3678794 


0. 0044048 


2. 0000000 


0.2968013 


0.1248231 


0.1353353 


0.0105122 


3.0000000 


0.2556222 


0.0349784 


0.0497871 


0.0148087 


4.0000000 


0.2401217 


0. 0011687 


0.0183156 


0.0171469 


5.0000000 


0.2342887 


-0.0115542 


0.0067379 


0.0182921 


6.0000000 


0.2320936 


-0.0163420 


0.0024788 


0. 0188207 


7. 0000000 


0.2312676 


-0.0181437 


0.0009119 


0.0190556 


8.0000000 


0.2309567 


-0.0188217 


0.0003355 


0.0191571 


9.0000000 


0.2308398 


-0. 0190768 


0.0001234 


0.0192002 


10. 0000000 


0.2307958 


-0.0191728 


0.0000454 


0.0192182 



* I NPUT DATA SET USED* 

H = 1.0000000 
ICON = 1 

**RE SUIT S OF ADAM2 PREDICTOR-CORRECTOR METHOD** 



X 


YP 


YC 


YEXACT 


ERROR 


0.0000000 






1.0000000 




1.0000000 


0.3750000 




0.3678794 


-0.0071206 


2.0000000 


0.3125000 


0.1015625 


0.1353353 


0.0337728 


3.0000000 


0.3437500 


-0.0312500 


0.0497871 


0.0810370 


4.0000000 


0.3281250 


-0.0996094 


0.0183156 


0.1179250 


5. 0000000 


0.3359375 


-0. 1328125 


0.0067379 


0.1395504 


6.0000000 


0.3320313 


~*0 . 1 4990 2 3 


0.0024788 


0. 1523811 


7.0000000 


0.3339844 


-0. 1582031 


0.0009119 


0.1591150 


8.0000000 


0.3330078 


-0.1624756 


0.0003355 


0.1628110 


9.0000000 


0.3334961 


-0.1645508 


0.0001234 


0.1646742 


10. 0000000 


0.3332520 


-0. 1656189 


0.0000454 


0.1656643 
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PROGRAM SEVEN 



THIRD ORDER ADAMS PREDICTOR-CORRECTOR METHOD. ODE: 
Y* =— Y WITH Y ( 0 ) =1 . TO SOLVE ANOTHER ODE CHANGE 
INITIAL CONDITIONS AND FUNCTION SUBROUTINES. 



J»> iJU OU* X V' X -X* X X X X X X X X X X X X X X X X X X X V* V*' X X X O/ X x x x x x x X X X X x X X# X X 
^ 'c ^ ^ '.'*■* *^ 'i ^ ^ ^ ^ v V r t ^ ^ ^ 'f ' t ^ ' -r 'r v -r* -r* T ^ -t* 



FORTRAN PROGRAM FOR THIRD ORDER ADAMS PREDICTOR- 
CORRECTOR METHOD IN THE NUMERICAL SOLUTION OF THE 
ORDINARY DIFFERENTIAL EQUATION IN THE FORM OF 
DY/DX=F(X,Y) WITH THE INITIAL CONDITION Y(X0)=Y0 
THE PARAMETERS TO THE PROGRAM HAVE THE FOLLOWING 
MEANINGS.. . 

XMAX _THE TERMINAL BOUNDARY CONDITION 

EPS THE CONVERGENCE TEST CONSTANT 

MAX THE MAXIMUM NUMBER OF ITERATIONS 

YP THE PREDICTED VALUE OF THE NEXT 

SOLUTION POINT 

YC _THE CORRECTED VALUE OF THE NEXT 

SOLUTION POINT 

FP THE FUNCTION EVALUATION 

H THE STEP SIZE TO BE USED TO ADVANCE 

THE POINT OF INTEGRATION 

ICON THE DESIRED POINT OF PRINTING THE 

COMPUTED SOLUTION FOR EVERY FIXED 
NUMBER OF INTEGRATION POINTS 
THE TRUE SOLUTION VALUES 



YEXACT 
ERROR , 

FCT. 

EXACT 



THE ABSOLUTE DEVIATION 
SOLUTION FROM THE TRUE 
SUBROUTINE TO 



FUNCTION 

FUNCTION 
FUNCTION 
SOLUTION 



OF THE NUMERICAL 
SOLUTION 
EVALUATE ODE 



TO COMPUTE TRUE 



SUBROUTI NE 
VALUES 

* RKUTTA THE SUBROUTINE USED TO GENERATE 

* NEEDED STARTING VALUES 






INITIALIZE INPUT CONSTANTS 

XMAX= I 0.0 
MAX= 2 

EPS=0. 000001 
K=2 

READ STEP SIZE AND CONTROL PRINTING 

5 READ ( 5 , 100 ) H,ICON 

SPECIFY INITIAL CONDITION 
X0=0.0 
Y0=1 . 0 
YR = Y0 

TEST FOR END OF INPUT DATA 

IF ( ICON— 0 ) 45,45,7 

PRINT INPUT DATA SET 

7 WRITE16, 1000) H,ICON 

WRITEI6 ,1001 ) XO , Y 0 

COMPUTE NEXT SOLUTION POINT 



193 



if- if- if- if- if- if- if- if if if if if if if if if if if if if if if if if if if if if if 



non non non non non ooo non 



NC=( XMAX+H/2. ) /H 
FO=FCT( X0,Y0) 

X 1 =XO+H 
X2=X0+2.*H 
X=X0+3 .*H 

CALL RKUTTA(X0,YR,Y1,Y2,K,H) 

F1=FCT(X1,Y1 ) 

F2=FCT1X2, Y2) 

DO 40 N=3,NC 
M= 1 

YP=Y1+H/12.*(23.*F2-16.*F1+5.*F0) 

FP=FCT(X,YP) 

11 YC=Y2+H/12.*( 5.*FP+8.*F2^-F1) 

DELY=ABS (YC-YP ) 

TEST FOR CONVERGENCE 

IF(DELY-EPS) 30,30,15 

TEST IF MAXIMUM NUMBER OF ITERATIONS IS SATISFIED 

15 I F ( M— MAX) 20,20,30 

20 FP=FCT(X,YC) 

M=M+ 1 
GO TO 11 

COMPUTE TRUE SOLUTION VALUE 
30 YE XACT = E XACT ( X ) 

COMPUTE ERROR IN CALCULATION 
ERROR= YEXACT-YC 

TEST IF DESIRED POINT OF PRINTING IS REACHED 

I F ( (N/ I CON* ICON) .EQ.N) WRI T E (6 , 100 2 ) X , YP , YC , YEXACT , 

1 ERROR 

UPDATE REQUIRED PARAMETER VALUES 

F0=F 1 
FI =F2 
F2 = FP 
Y 0= Y 1 
Y1=Y 2 
Y2=YC 
X= X+H 

40 CONTINUE 

LOOP BACK TO READ NEXT SET OF INPUT DATA 

GO TO 5 
45 STOP 

100 FORMAT (FLO. 7, I 10) 

1000 FORMAT (////41X, • *INPUT DATA SET US ED* ' //43X, ' H =', 
1F10.7/43X, 'ICON = • ,110) 

1001 FORMAT (//29X,*** RESULTS OF ADAM3 PREDICTOR-CORRECTOR 

1 METHOD** 1 // 10X , ' X ',5X,' YP ’,5X,' YC ', 

15X , ' YEXACT ' , 5X , • ERROR • // 12 X , F 1 0 . 7 ,44X , 

1 FI 0 . 7 ) 

1002 FORMAT ( '0',9X,F15.7,2X,F15.7,2X,F15.7,2X,F15.7,2X, 

IF 15 » 7 ) 

END 



SUBROUTINE RKUTTAI XX,YY,Y1 ,Y2, KK,HH) 
NN = 1 

17 CK1=HH*FCT(XX,YY) 

CK2=HH*FCT( XX+HH/2. 0,YY+CK 1/2.0) 
CK3=HH*FCT (XX+HH/2.0, YY+CK2/2.0 ) 
CK4=HH*FCT { XX+HH , Y Y+CK3 ) 
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GO TO ( 101,102} ,NN 

101 Yl=YY+(CKl+2.0*CK2+2.0*CK3+CK4)/6.0 
YY=Y1 

77 XX=XX+HH 

YEXACT=EXACT (XX ) 

PRRHC-YP VAT T— V Y 

WRITE16, 1004) XX,YY,YEXACT , ERROR 
NN=NN+1 

IF { NN-KK } 17,17,99 

102 Y2=YY+(CK1+2.0*CK2+2.0*CK3+CK4)/6.0 
YY=Y2 

GO TO 77 

1004 FORMAT (»0',9X,F15.7,2X,F15.7,16X,F15.7,2X,F15.7) 
99 RETURN 
END 



FUNCTION FCT ( XN , YN ) 

FCT=- YN 

RETURN 

END 



FUNCTION EXACT(XK) 

EXACT=EXP(-XK) 

RETURN 

END 
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♦INPUT DATA SET USED* 

H = 0.5000000 
ICON = 2 

♦♦RESULTS OF ADAM3 PREDICTOR-CORRECTOR METHOD** 



X 


YP 


YC 


YEXACT 


ERROR 


0.0000000 






1.0000000 




0.5000000 


0.6067709 




0.6065306 


-0.0002403 


1.0000000 


0.3681710 




0.3678794 


-0.0002916 


2.0000000 


0.2630881 


0. 1307259 


0.1353353 


0.0046093 


3.0000000 


0.0949250 


0.0457278 


0.0497871 


0.0040593 


4.0000000 


0.0333163 


0.0160214 


0.0183156 


0.0022942 


6.0000000 


0.0041118 


0.0019661 


0.0024788 


0.0005126 


7. 0000000 


0.0014410 


0. 0006887 


0.0009119 


0.0002232 


8.0000000 


0.0005048 


0.0002412 


0.0003355 


0.0000942 


9.0000000 


0.0001768 


0.0000845 


0.0001234 


0. 0000389 


10.0000000 


0. 0000619 


0.0000296 


0.0000454 


0.0000158 



♦INPUT DATA SET USED* 

H = 1.0000000 
ICON = 1 

♦♦RESULTS Or ADAM3 PREDICTOR-CORRECTOR METHOD** 



X 


YP 


YC 


YEXACT 


ERROR 


0.0000000 






1.0000000 




1.0000000 


0.3750000 




0.3678794 


-0.0071206 


2.0000000 


0.1406250 




0.1353353 


-0. 0052897 


3.0000000 


0.1888022 


0.0454788 


0.0497871 


0.0043082 


4.0000000 


0.0217022 


0.0021872 


0.0183156 


0.0161284 


5.0000000 


0.0785821 


-0. 0024490 


0.0067379 


0.0091869 


6.0000000 


-0.0525024 


-0.0057783 


0.0024788 


0.0082570 


7. 0000000 


-0.0479046 


0. 0015025 


0.0009119 


-0. 0005906 


8. 0000000 


-0.0577729 


-0.0018528 


0.0003355 


0.0021883 


9.0000000 


0.0527026 


0.0029584 


0.0001234 


-0.0028350 


10.0000000 


^0.0540226 


-0. 0020290 


0.0000454 


0.0020744 
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PROGRAM EIGHT 



FOURTH ORDER ADAMS PREDICTOR-CORRECTOR METHOD. ODE 
Y' =- Y WITH Y( 0) =1. TO SOLVE ANOTHER ODE CHANGE 
INITIAL CONDITIONS AND FUNCTION SUBROUTINES. 



C 

c 

c 

c 

c 

c 



c 

c 

c 

c 



c 

c 

c 

t 

c' 

c 

c 



c 

c 

c 

c 

c 

c 



c 

c 



i£5(r # sje # £ st Jf: # $ £ # # % sjc # £ * # ^ & # # # # =Jc jf: £ # ^ ^ ^ s^c # * # ^ ajc 3 

FORTRAN PROGRAM FOR FOURTH ORDER ADAMS PREDICTOR- 
CORRECTOR METHOD IN THE NUMERICAL SOLUTION OF THE 
ORDINARY DIFFERENTIAL EQUATION IN THE FORM OF 
DY/DX=F(X,Y) WITH THE INITIAL CONDITION Y(XO)=YO 
THE PARAMETERS TO THE PROGRAM HAVE THE FOLLOWING 

THE TERMINAL BOUNDARY CONDITION 
THE CONVERGENCE TEST CONSTANT 
THE MAXIMUM NUMBER OF ITERATIONS 
THE PREDICTED VALUE OF THE NEXT 
SOLUTION POINT 

THE CORRECTED VALUE OF THE NEXT 
SOLUTION POINT 
THE FUNCTION EVALUATION 
THE STEP SIZE TO BE USED TO ADVANCE 
THE POINT OF INTEGRATION 
THE DESIRED POINT OF PRINTING THE 
COMPUTED SOLUTION FOR EVERY FIXED 
NUMBER OF INTEGRATION POINTS 
TRUE SOLUTION VALUES 






c 


* 


MEANINGS 


m m m 


c 


❖ 


XMAX 




c 


* 


EPS 




c 


❖ 


MAX 




c 


* 


YP 




c 

c 


* 


YC, 




c 


❖ 






c. 




FP 




c 


* 


H 




c 


* 






c 


* 


ICON 




c 






1 


c 


•JU 

•r- 






c 


* 


YEXACT 




c 


* 


ERROR 




c 








c 


* 


FCT 




c 


*r 






c 


* 


EXACT 




c 


* 






c 


❖ 


RKUTTA. 




c 






'1 


c 









OF THE NUMERICAL 
SOLUTION 
EVALUATE ODE 



THE ABSOLUTE DEVIATION 
SOLUTION FROM THE TRUE 
FUNCTION SUBROUTINE TO 
FUNCTION 

FUNCTION SUBROUTINE TO COMPUTE TRUE 
SOLUTION VALUES 

THE SUBROUTINE USED TO GENERATE 
NEEDED STARTING VALUES 

#*£****# S *** j*#**#* i* *:}:**:(£*#** ***:£* ** # 



INITIALIZE INPUT CONSTANTS 

XMAX= 10.0 
MAX=2 

EPS=0. 000001 
K=3 

READ STEP SIZE AND CONTROL PRINTING 

READ (5,100) H» I CON 

SPECIFY INITIAL CONDITION 

X0=0.0 
Y0=1 . 0 
YR=Y0 

TEST FOR END OF INPUT DATA 

I F ( I CON-O ) 45,45,7 

PRINT INPUT DATA SET 

WRITE! 6,1000) H,ICON 
WRITE(6, 1001) XO , Y 0 

COMPUTE NEXT SOLUTION POINT 
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ik * # X -H- * # # -M- H 



non non non non non non ooo 



NC=(XMAX+H/2. )/H 
FO=FCT(XO,YO) 

X1=X0+H 
X2=X0+2 • *H 
X3=X0+3.*H 
X=X0+4 «*H 

CALL RKUTTA(X0,YR,Y1,Y2,Y3,K,H) 

F1=FCT( XI, Yl) 

F2=FCT (X2 t Y2) 

F3=FCT (X3 , Y3 ) 

DO 40 N=4 r NC 
M= 1 

YP=Y3+H/24.*(55.*F3-59.*F2+37.*F1-9.*F0) 

FP=FCT(X,YP) 

11 YC=Y3+H/24.*(9.*FP+19.*F3-5.*F2+F1 ) 

DELY=ABS(YC-YP ) 

TEST FOR CONVERGENCE 

I F ( DELY-EPS) 30,30,15 

TEST IF MAXIMUM NUMBER OF ITERATIONS IS SATISFIED 

15 I F ( M— MAX ) 20,20,30 

20 FP=FCT (X , YC) 

M=M+ 1 
GO TO 11 

COMPUTE TRUE SOLUTION VALUE 
30 Y EX ACT = EXACT ( X ) 

COMPUTE ERROR IN CALCULATION 
ERR0R= YEXACT-Y C 

TEST IF DESIRED POINT OF PRINTING IS REACHED 

IF ( (N/IC0N*IC0N).EQ.N) WRI T E (6 , 1 00 2 ) X ,YP, YC, YEXACT , 
1 ERROR 

UPDATE REQUIRED PARAMETER VALUES 

Y0=Y 1 
Yl =Y 2 
Y2=Y3 
Y3=YC 
F0 = F 1 
F1 = F2 
F2=F3 
F3 = FP 
X= X+H 

40 CONTINUE 

LOOP BACK TO READ NEXT SET OF INPUT DATA 

GO TO 5 
45 STOP 

100 FORMAT (F10.7 , 1 10) 

1000 FORMAT (////41X, '*INPUT DATA SET USED* ' //43X , 'H =' 

1F1G. 7/43X, • ICON =' ,110) 

1001 FORMAT ( / /29X, ' * s * RESULTS OF ADAM4 PREDICTOR-CORRECTOR 
1METHCD** 1 // 10X, • X '^X,* YP '^X,* YC ' 

1 5X , 1 YEXACT ' , 5X , • ERROR • // 12X , FI 0 . 7 , 44X , 
1F10.7) 

1002 FORMAT ( • 0 • ,9X , FI 5- 7 , 2X , F 15 .7 , 2X , FI 5 . 7, 2X , F 15 . 7 , 2X , 

IF 1 5. 7 ) 
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SUBROUTINE RKUTTAI XX , YY ,Y1 ,Y2 , Y3 ,KK t HH) 

NN=1 

CK1=HH*FCT(XX,YY) 

CK2= HH^FCT ( XX+HH /2 . 0 , YY+CK 1 / 2 . 0 ) 

CK3=HH*FCT ( XX+HH /2.0, YY+CK 2/2.0) 

CK4=HH*FCT (XX+HH, YY+CK3) 

GO TO ( 101,102,103) , NN 

101 Yl=YY+(CKl+2.0*CK2+2.0*CK3+CK4)/6.0 
YY=Y 1 

77 XX=XX+HH 

YEXACT=EXACT (XX) 

ERROR=YE XACT— YY 

WRITEI6, 1004) XX, YY,YEXACT, ERROR 
NN=NN+ 1 

IF(NN-KK) 17,17,99 

102 Y2=YY+(CKl+2.0*CK2+2.0*CK3+CK4)/6. 0 
YY=Y2 

GO TO 77 

103 Y3- YY+ ( CK 1+2 . 0*CK2+2 . 0 *CK3+CK4 ) /6. 0 
YY=Y3 

GO TO 77 

1004 FORMAT ( • 0 • ,9X , F 15. 7 , 2X , F15 . 7 , 16X,F 15.7 ,2X, F 1 5. 7) 
99 RETURN 
END 



FUNCTION FCT (XN, YN ) 

FCT=-YN 

RETURN 

END 



FUNCTION EXACT(XK) 
EXACT = EXP (-XK) 
RETURN 
END 
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* I N PUT DATA SET USED* 

H = 0.5000000 
ICON = 2 

**RE SULTS OF ADAM4 PREDICTOR-CORRECTOR METHOD** 



X 

0.0000000 


YP 


YC 


YEXACT 

1.0000000 


ERROR 


0.5000000 


0.6067709 




0 .6065306 


-0.0002403 


1.0000000 


0.3681710 




0.3678794 


-0.0002916 


1. 5000000 


0.2233955 




0.2231301 


-0.0002654 


2.0000000 


0.1397457 


0.1352788 


0.1353353 


0.0000565 


3.0000000 


0.0512678 


0. 0495746 


0.0497871 


0.0002124 


4. 0000000 


0.0187947 


0.0181669 


0 .0183156 


0.0001488 


5.0000000 


0 . C068879 


0.0066569 


0.0067379 


0. 0000810 


6.0000000 


0.0025232 


0.0024393 


0.0024788 


0.0000394 


7.0000000 


0.0009245 


0.0008938 


0.0009119 


0.0000180 


8.0000000 


0.0003387 


0. 0003275 


0.0003355 


0.0000079 


9.0000000 


0.0001241 


0.0001200 


0.0001234 


0.0000034 


10.0000000 


0.0000455 


0.0000440 


0.0000454 


0.0000014 





* INPUT 


DATA SET 


USED* 






H 


= 1.0000000 








ICON 


= 


1 






**RESULTS OF ADAM4 


PREDICTOR- 


•CORRECTOR METHOD** 


X 


YP 


YC 




YEXACT 


ERROR 


0.0000000 








1 .0000000 




1.0000000 


0.3750000 






0.3678794 


-0.0071206 


2.0000000 


C. 1406250 






0.1353353 


-0.0052897 


3.0000000 


0.0527344 






0.0497871 


-0.0029473 


4. 0000000 


0.0744628 


0. 0149522 




0.0183156 


0.0033634 


5.0000000 


0.C091045 - 


0.0007950 




0.0067379 


0.0075330 


6.0000000 


0.0319238 - 


0.0004662 




0.0024788 


0.0029450 


7.0000000 


-0.0306429 - 


■0.0027268 




0.0009119 


0.0036387 


8.0000000 


0.0368799 


0.0015700 




0.0003355 


-0.0012345 


9.0000000 


-0.0442679 - 


0.0027738 




0.0001234 


0.0028972 


10.0000000 


0.0550162 


0.0028112 




0.0000454 


-0.0027658 
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